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We present the implementation of a fast estimator for the full dark matter bispectrum of a three- 
dimensional particle distribution relying on a separable modal expansion of the bispectrum. The 
computational cost of accurate bispectrum estimation is negligible relative to simulation evolution, 
so the isotropic bispectrum can be used as a standard diagnostic whenever the power spectrum 
is evaluated. As an application we measure the evolution of gravitational and primordial dark 
matter bispectra in A/'-body simulations with Gaussian and non-Gaussian initial conditions of the 
local, equilateral, orthogonal and flattened shape. The results are compared to theoretical models 
using a 3D visualisation, 3D shape correlations and the cumulative bispectrum signal-to-noise, all 
of which can be evaluated extremely quickly Our measured bispectra are determined by (9(50) 
coefficients, which can be used as fitting formulae in the nonlinear regime and for non-Gaussian initial 
conditions. In the nonlinear regime with k < 2/iMpc~^, we find an excellent correlation between 
the measured dark matter bispectrum and a simple model based on a 'constant' bispectrum plus 
a (nonlinear) tree-level gravitational bispectrum. In the same range for non-Gaussian simulations, 
we find an excellent correlation between the measured additional bispectrum and a constant model 
plus a (nonlinear) tree-level primordial bispectrum. We demonstrate that the constant contribution 
to the non-Gaussian bispectrum can be understood as a time-shift of the constant mode in the 
gravitational bispectrum, which is motivated by the one-halo model. The final amplitude of this 
extra non-Gaussian constant contribution is directly related to the initial amplitude of the constant 
mode in the primordial bispectrum. We also comment on the effects of regular grid and glass initial 
conditions on the bispectrum. 
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I. INTRODUCTION 

Observations of the cosmic microwave background 
(CMB) [1-4] and large scale structure (LSS) [5, 6] are 
currently consistent with Gaussian primordial cosmolog- 
ical perturbations. Much effort has been undertaken in 
recent times to constrain a possible non-Gaussian con- 
tribution. Detection of a significant primordial bispec- 
trum or trispectrum would have major implications for 
the mechanism of inflation, possibly ruling out the sim- 
plest paradigm of canonical single-field slow-roll infla- 
tion. Given this possibility, the importance of developing 
methods that discriminate between the plethora of infla- 
tionary models is clear. Such general methods have been 
developed in the case of the CMB in [2, 4, 7-9]. This 
work exploited the use of a separable expansion of the 
underlying bispectrum or trispectrum in order to greatly 
reduce the computational cost involved in analysing gen- 
eral shapes. 

The CMB bispectrum is currently the most power- 
ful and direct probe of primordial non-Gaussianity, and 
higher resolution temperature and polarisation data will 
soon be available from Planck. Nevertheless, interest in 
the possibility of using observables of large scale struc- 
ture to test primordial non-Gaussianity has undergone 
a recent resurgence, due in part to the potentially three 
dimensional nature of the dataset as a result of its red- 
shift dependence. In particular, the scale-dependent bias 
induced by non- Gaussian local initial conditions on the 
halo power spectrum has been shown in [10] to offer an 
additional and powerful probe. Much work has been un- 
dertaken to improve analytic predictions for the impact of 
non- Gaussian initial conditions on the matter and galaxy 
power spectrum and bispectrum (see for example [11-21] 
and references therein). In [6, 22] constraints on local- 
type non-Gaussianity were found at levels competitive 
with those obtained using CMB data. Despite these ad- 
vances it is notable that until recently only local type 
non-Gaussianity had been studied using large-scale struc- 
ture, owing to the difficulty in generating generic initial 
conditions. This situation greatly improved due to work 
by Wagner and Verde [23, 24] and work by Scoccimarro 
et al. [25]. However, both of these approaches become in- 
efficient for non-separable bispectra. In [26] the possibil- 
ity of using the separable expansion method, exploited to 
great effect in the case of the CMB [2], was explored. This 
approach was verified in [27] to allow for a far more effi- 
cient generation of non-Gaussian initial conditions than 
had previously been possible in the literature. Further- 



more, the separable decomposition methodology allows 
for the study of non-separable shapes. This has allowed 
for a robust and truly general approach towards the gen- 
eration of primordial non-Gaussian initial conditions for 
use in AT-body simulations. In this paper we exploit this 
approach to set up and run TV-body simulations with 
non-Gaussian initial conditions. In particular, we study 
initial conditions given by local, equilateral, orthogonal 
and flattened bispectra, respectively. We choose the flat- 
tened (or trans-Planckian) model as an explicit example 
of an inherently non-separable bispectrum. We exploit 
the modal method to efficiently and accurately recon- 
struct the full 3D matter bispectrum at each redshift of 
interest for each of these non- Gaussian models. This al- 
lows us to correlate the results of the N-hodj simulations 
with analytic estimates, allowing us to identify clearly the 
regime at which nonlinear corrections become important. 
This applies first to accurately determining the gravita- 
tional bispectrum well into the nonlinear regime, which is 
necessary in order to differentiate the more subtle impact 
of primordial non-Gaussianity. The correlation measure 
proves useful in distinguishing the different primordial 
shapes. We also thoroughly test the impact of start- 
ing TV-body simulations using either glass or regular grid 
initial conditions. We emphasise that our purpose here 
is to study the detailed nature of the underlying mat- 
ter bispectrum derived from A/'-body simulations. These 
methods can be equally applied to efficiently extract the 
galaxy bispectrum from huge survey data sets or to pre- 
dict the halo bispectrum from simulations in a realistic 
observational context, but this will be the subject of fu- 
ture work. 

The paper is organised as follows. In Section II we 
review the generation of non-Gaussianity due to nonlin- 
ear gravitational evolution, together with the physically- 
motivated models we will use to describe our results. In 
Section III we review primordial non-Gaussianity, intro- 
duce the specific models studied in this paper and dis- 
cuss their impact on the matter bispectrum, proposing 
a new time-shift model. In Section IV we discuss the 
bispectrum estimation methodology based on a separa- 
ble expansion which allows extremely efficient estimation 
of the full A'-body bispectrum, as well as the cumula- 
tive signal-to-noise and its correlation to theoretically- 
predicted bispectra. We discuss the simulation setup, 
impact of glass initial conditions, validations and conver- 
gence tests in Section V. In Section VI we present our 
results for the gravitational bispectrum and the simple 
fits thereof, while primordial bispectra in non- Gaussian 
simulations and their fits are discussed in section VII. Fi- 
nally, in Section VIII we present our conclusions. Readers 
familiar with the modal methodology and interested pri- 
marily in the measured bispectra and testing of various 
fitting formulae may wish to start with Section VI and 
follow the references to the earlier sections given there. 
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II. THE DISTRIBUTION OF MATTER 
Power spectrum and transfer functions 

The distribution of matter in the universe can be de- 
scribed by its fractional overdensity ^(x) = (p(x) — 
where p is the spatial average of the density p(x). An 
important prediction of cosmological models is the prob- 
ability of finding a certain configuration (5(x) in the uni- 
verse. The simplest possibility is a Gaussian pdf, which 
is determined by the two-point correlation function, or 
in Fourier space by the power spectrum P^, 

(J(ki)(^(k2)) = (2^)'(5z,(ki + k2)P,(^i). (1) 

Here we assumed statistical homogeneity, leading to the 
Dirac delta function, and statistical isotropy, which im- 
plies that the power spectrum only depends on the mag- 
nitude of the wavevector. 

Perturbations from infiation are usually described by 
the comoving curvature perturbation 7^, which is related 
to the primordial potential $ by <l> = —37^/5, in lin- 
ear perturbation theory during matter domination. The 
density perturbation 5 can be obtained using the linear 
Poisson equation 

^(k,.) = I ^^^^^^ *(k) ^ M(fc,^)$(k), (2) 

where T{k) is the linear transfer function at low redshift 
normalised to T{k) = 1 on large scales and calculable 
with CAMB [28], and D{z) is the linear growth function 
for r^rad = [29] normalised to D{z) = 1/(1 + z) during 
matter domination. Later we will also find it convenient 
to use the growth function normalised to unity today, 
D{z) = D{z)/D{Q). 

At late times and on small scales, the linear treatment 
breaks down and A^-body simulations must be used to 
find the nonlinear transfer of primordial to late time per- 
turbations. A widely used fit of the nonlinear power spec- 
trum is HALOFIT [30], which is included in CAMB [28]. 
However, in our own bispectrum analysis we will gener- 
ally use the actual nonlinear power spectrum measured 
from the A'-body simulations. 

Matter bispectrum 

At late times, the model of a Gaussian pdf for the 
density perturbation is over-simplified because nonlinear 
gravitional collapse will lead to higher order correlations. 
To study deviations from Gaussianity, it is useful to con- 
sider the bispectrum B§ (or three-point correlator trans- 
form), which is defined by 

{S{ki)S{k2)S{ks)) = {27rf Sd (S.k,) Bsih, ^2, fe), (3) 



assuming again statistical homogeneity and isotropy. 
The Dirac delta function imposes a triangle constraint 
on the bispectrum arguments /ci, /c2, /ca. The space of tri- 
angle configurations is sometimes called a tetrapyd [8] 
and is shown in Fig. 1. Additionally to the triangle 
constraint, we impose a cut off at /cmax, corresponding 
to the smallest scale under consideration. The bispec- 
trum shape, i.e. its dependence on the wavenumbers /c^, 
gives important information about the physics which in- 
duced the bispectrum. Not only does it help to disentan- 
gle late time non- Gaussianity, e.g. induced by nonlinear 
gravitational collapse, from primordial infiationary non- 
Gaussianity, but it also opens the intriguing possibility of 
distinguishing between different inflationary models. To 
see this explictly we will review the perturbative treat- 
ment of such non-Gaussianities in the following sections. 



Tree- level gravitational matter bispectrum 

The equations of motion for a perfect pressureless fluid 
in a homogeneous and isotropic universe contain nonlin- 
earities, which induce a non-vanishing bispectrum (see 
e.g. [31, 32] for details). Speciflcally these are the V-(^v) 
term in the continuity equation and the (v • V)v term 
in the Euler equation, where v is the peculiar velocity, 
which describes deviations from the background Hubble 
flow. To see how these terms induce a non-Gaussian den- 
sity perturbation one makes a perturbative ansatz for the 
density perturbation (see [31] for a review) 

CO 

d{k,t) = J2a{trSn{k) (4) 




Figure 1. Space of triangles with sides A;i,A:2,/c3, i.e. each 
point inside the tetrapyd volume corresponds to a triangle 
configuration. Squeezed, folded and equilateral configurations 
are highlighted. 
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and similarly for the peculiar velocity, which is assumed 
to be described completely by its divergence, = V - 
V. Then the equations of motion imply (see [31] and 
references therein) 

(5n(k) = j d^qi - j c^^qnfo(k - qi qn) 

i^^'Hqi,...,qn)^i(qi)---^i(qn), (5) 

is) 

and determine the kernels Fk°\ e.g. 

F,^Hk„k,) = - + -— + + j . 

(6) 

Here the first and second terms come, respectively, from 
(5 V • V and v • V5 in the continuity equation, while the 
last term comes from (v • V)v in Euler's equation [32]. 

At high redshift, the expansion (4) is dominated by 
5i and 52 - F^'hl. If we assume ^i to be Gaussian, 
the leading order contribution to the gravitational bis- 
pectrum is given by [33] 

Br''{kuk2, h) = 2Pt{h)Pt{k2)F^'\\iiM) + 2 perms, 

(7) 

where is the power spectrum of the linear perturba- 
tion 6i. 

To discuss the shape of this gravitational bispec- 
trum let us consider two-dimensional slices through the 
tetrapyd shown in Fig. 1, with ki -\- k2 -\- k-^ = const. We 
denote one side of these two-dimensional triangular slices 
as a and parameterise the other direction with A 
slice through the gravitational bispectrum (7) is shown 
in Fig. 2a (see [34] for more slices at slightly different 
length scales). The bispectrum is maximal at the edges 
of the plot, corresponding to flattened triangle config- 
urations, where /ci + /c2 = or permutations thereof, 
i.e. where the wavevectors are parallel or anti-parallel 
to each other. However there is a suppression in the cor- 
ners of the plot, corresponding to squeezed triangle con- 
figurations with ki <^ k2 ^ ks or permutations thereof. 
Non-flattened and equilateral triangle configurations in 
the centre of the plot are also suppressed. This tree level 
gravitational bispectrum is also illustrated in Fig. 3 as a 
function over the full tetrapyd. 

To understand the basic shape of the gravitational 
bispectrum (7) we plot in Fig. 2b the expression 
2P^{ki)P^{k2) + 2 perms, which corresponds to replac- 
ing ^2*^ in (7) by a constant. Comparing Fig. 2b with 
Fig. 2a shows that the configuration dependence of the 



Details can be found in [7]. The slice parameters a, /3 should not 
be confused with the expansion coefficients aj^'^^ and 
to be defined later. 



kernel (6), which is induced by the terms containing 
scalar products, leads to an enhancement of flattened 
and particularly folded configurations, where two of the 
wavevectors k^ equal each other. Non- flattened configu- 
rations are relatively suppressed. As we go along the edge 
of the plot in Fig. 2a, from folded {ki = k2 = ks/2 or per- 
mutations) to squeezed configurations {ki <C /c2 ~ /cs or 
permutations), the bispectrum shape reflects the shape of 
the power spectrum P/', which peaks at /Ceq ~ 0.02/i/Mpc 
and then decreases with decreasing ki because P^{ki) oc 
ki on large scales. 

Further discussion is required in the squeezed limit, 
where k2 ~ — ks. Let us consider the regime where 
ki < /Ceq and /c2, > /Cgq. The term with P^{k2)P^{k^) 
in (7) is small since the small scale power spectra decrease 
rapidly with increasing /c2,/c3 and F^^^ {)<.2^^?) vanishes 
for k2 = — ks. The other two permutations in (7) de- 
pend on the angle between the large-scale wavevector 
ki and the small-scale wavevectors k2,k3. First con- 
sider the case where the large scale is approximately 
perpendicular to the two small scales, i.e. ki • k2 ~ 
-ki • k3 ^ 0. Then Fi^\ki,k2) ^ Fi'^(ki,k3) ^ 5/7 
and P^{ki)P^{k2) and P^{ki)P^{k^) decrease as we ap- 
proach more squeezed triangles, implying a suppressed 
bispectrum in the squeezed limit. In the other limit, 
the large-scale wavevector ki is not orthogonal to the 
small scale wavevectors, and so is aligned with one or 
other of k2, k3. Then the squeezed limit /ci ^ implies 
F2^'^(ki,k2) OC k:[^ oo and F2^'^(ki,k3) ex -k^^ 
— oo, which can be seen in Fig. 2c. However, in the sum 

{k^)Pt {k2)F^'^ (ki , k2) + P/- (fci )P/' {k3)F^'^ (ki , ks) , 

the two terms containing k^^ divergences in the kernel 
approximately cancel, because k2 ~ ks and ki • k2 ~ 
— ki • k3. Also in the limit /ci ^ 0, the divergences of the 
kernels are regulated by the large scale power spectrum 
because k^^P^{ki) = const, on very large scales. Fig. 2a 
shows that the divergences are indeed cancelled in the 
total bispectrum (7) and the squeezed limit is suppressed 
if the large-scale wavenumber satisfies ki < /Ceq. 



Gravitational matter bispectrum beyond tree level 

Loop corrections 

The tree level prediction for the gravitational matter 
bispectrum (7) is only a good approximation on large 
scales and can be improved by including so-called loop 
corrections, which were derived for Gaussian initial con- 
ditions in [32] and extended to include non- Gaussian ini- 
tial conditions in [11]. Important loop corrections can be 
included simply by replacing the linear power spectrum 
P^ by the nonlinear power spectrum Ps in the tree level 
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(a) (b) (c) 

Figure 2. Two-dimensional slices of (a) the leading order gravitational bispectrum from (7), (b) the expression 2P^ (/ci)P^ (^2) + 

perms and (c) the kernel ^2*^* defined in (6) (with cut-offs in the squeezed limit where the kernel diverges). All slices are at fixed 
ki -\- k2 + ks = 0.24/i/Mpc and we restrict the plot to ki > 0.0013/i/Mpc, corresponding to a 5000Mpc//i box. The linear power spectra 
are evaluated at 2; = 30. 



expression (7) for Gaussian initial conditions, that is, 
Bl^^l = 2Ps{h)Ps{k2)F^'\k,,k2) + 2 perms . (8) 

We will find the expression (8) to be a key approxima- 
tion to the full gravitational bispectrum as we probe be- 
yond the mildly nonlinear regime later. Of course, this 

(s) 

omits several loop corrections containing powers of F2 
or higher order kernels, but it is easy to evaluate with the 
nonlinear power spectrum from C AMB [28] . We will dis- 
cuss further loop corrections much more quantitatively 
in a forthcoming publication [35]. 



In detail, the halo model prediction for the dark matter 
bispectrum is given by three contributions corresponding 
to the three cases that the three points of the three point 
function lie in one, two or three halos (see e.g. [36-38]): 

^HM = ^IH + ^2H + ^3H, (9) 

where 

^ih(^i,^2,^3) = ^ / dmn{m)Y[p{ki,m), (10) 



Halo models 

At sufficiently small scales, well probed by our simula- 
tions, the perturbative treatment breaks down and sim- 
ulations and phenomenological models must be used. In 
the strongly nonlinear regime the halo model can be used 
as a phenomenological model for the dark matter distri- 
bution (see [36] for a review). Recently in Ref. [37] it was 
demonstrated that combining 1-loop perturbation theory 
with the halo model describes the matter bispectrum in 
simulations at the (9(10%) level on nonlinear scales for 
equilateral and isosceles bispectrum slices. Similar re- 
sults were obtained in Ref. [38], where local non-Gaussian 
initial conditions were also considered and compared with 
simulations for k < 0.3/i/Mpc at z < 1. In the mildly 
nonlinear regime, when the approximation that all dark 
matter particles are inside halos is no longer valid, the 
halo model becomes less accurate than in the strongly 
nonlinear regime. It should also be noted that important 
ingredients of the halo model like the halo density profile 
and halo mass function cannot be derived analytically 
but are obtained from fits to A/'-body simulations. 



^2H 



(/ci, /c2, /cs) = ^ J dmin{mi)p{ki,mi) J dm2 

n{m2)p{k2,m2)p{ks, m2)P^(/ci, mi, 7712) + perms, 

(11) 



B3u{ki,k2,k3) 



■ 3 

Y[ / dmin{mi)p{ki,mi) 



Bh{ki,mi; /c2, 1712; ks, ma). (12) 



Here n(m) is the mass function and p(/c, m) is the Fourier 
transform of the density profile of a halo of mass m, 
i.e. for spherically symmetric density profiles (like the 
commonly used NFW profile [39]) 

X A f 1 2 ^ .sm(kr) 
p[k^m) = 47r / drr pyr^m) — ^ . (13) 

Ph and Bh denote the halo power spectrum and bispec- 
trum, which can be related to the tree level dark matter 
power spectrum and bispectrum on large scales using bias 
relations. 
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Constant 'halo ' model 

As we consider the more strongly nonlinear regime, 
k > 1/i/Mpc at z = 0, the halo model bispectrum is 
dominated by the 1-halo contribution Bm. Neglecting 
the overall bispectrum amplitude, we find that in this 
regime the 1-halo shape for Gaussian initial conditions 
can be approximated by a simple constant bispectrum, 

^fZnst(^l' ^2, fc3)U = Ci D^-iz) ih +k2 + k^f . (14) 

Here, 'constant' refers to constancy on /ci + /c2 + /cs = 
const, slices (see [7]), noting that this is a bispectrum 
typical of isolated point-like objects (e.g. particles ran- 
domly distributed in a box have a completely constant 
shot noise bispectrum). The expression (14) has an over- 
all wavenumber scale-dependence with exponent v and a 
time-dependence on the linear growth function D{z) with 
halo exponent rih- The wavenumber scaling v ^ —1.7 is 
chosen such that it approximately refiects the scaling for 
equilateral triangle configurations in this regime as mea- 
sured in our simulations (and which is approximately pre- 
dicted by the halo model for Gaussian initial conditions 
[37, 38]). The exponent rih is similarly defined by the 
fast growth factor appropriate for the halo model for the 
scales under study 0.2/i/Mpc ^ k < 2/i/Mpc, typically 
with Uh ~ 6-8. For /cmax = 2/i/Mpc and z = 0, this 
simple model achieves a shape correlation of 99.7% with 
the 1-halo contribution Bm and 99.3% with the full halo 
model bispectrum (9). 

In later sections, we will note that (14) provides an ex- 
cellent approximation to the late-time bispectrum in the 
nonlinear regime when we use it in a simple fitting for- 
mula together with the modified tree-level gravitational 
bispectrum (8). We shall investigate the other halo con- 
tributions in more detail elsewhere [35]. 

As we shall see, on smaller nonlinear scales (large 
k > 1/i/Mpc) the fast bispectrum growth begins to slow 
down by the present day z = 0. In this case, we should 
really replace the power law growth D^^ (z) for the con- 
stant mode using a more general growth factor T{k, z, Zi), 
where the 'slice' or 'average' wavenumber 



^= (A:i-hA:2 + A:3)/3. 



(15) 



For future reference, it is convenient to use this growth 
rate in a more general integral form of the 'constant' 
model (14), that is, 

Bi:::nstikuk2, ^3)1. = nk z, zo B^uk, (w) 



f 



g{k,z)dzBl^l,,{k,zO, 



where Zi is the redshift at which this rapid 'halo' growth 
takes hold (it has an implicit k dependence). The quan- 
tity 5const(^5 ^i) represents the initial condition at z = Zi 
for the constant part of the bispectrum. Naively, we 



might take this to be B'^J^^,{k,zO = Bp^k,k,k)U=,,, 
that is, the equilateral or constant part of the tree-level 
gravitational bispectrum (7) at z = Zi. However, to 
date, determining the amplitude of this initial constant 
bispectrum has relied on simulations. We will use this 
simple model to characterise both the time- and scale- 
dependence of the gravitational bispectrum, as well as 
primordial non-Gaussianity as we approach the strongly 
nonlinear regime. 



Alternative phenomenological fit 

An alternative description of the gravitational matter 
bispectrum in the non-perturbative regime was proposed 
in Ref. [13], who constructed a fitting formula which in- 
terpolates between the perturbative prediction on large 
scales and a local-type bispectrum on small scales, which 
was suggested by early simulations. In their terminol- 
ogy, the small-scale bispectrum was denoted as a 'con- 
stant reduced' bispectrum, which implies that it takes the 
same form as the local shape Bs ~ Ps{ki)Ps{k2) -\- perms ^ 
in contrast to the constant shape (14) above. Recently 
Ref. [14] extended this fitting formula with updated sim- 
ulations into mildly nonlinear scales 0.03/i/Mpc < k < 
0.4/i/Mpc at < z < 1.5. 

We review for the sake of completeness this phe- 
nomenological 9-parameter fit, for which we will test the 
regime of validity. The linear power spectrum in (7) is 
replaced by the nonlinear one P^, and the symmetrised 

(s) 

kernel F2 is replaced by 



F2^^^^^(ki,k2) = -a(ni, /ci)a(n2,/c2) 



7 

+ ^ + ^ 1 b{ni,ki)b{n2, k2) 



2 kik2 \k2 ki 

2 / ki-k2 ^ ' 

7 V kik2 



c{ni,ki)c{n2,k2) , (17) 



where 
a(n, k) 
b{n, k) 



l+ag^(z)V0.7Q3(n)(aig)"^+" 
1 + (aig)''2+« 

1 + 0.2a3(n + 3)(ga7)"+^+''" 
1 + (ga7)3-5+"+«8 



(18) 
(19) 



_ l+4.5a4/[l-5 + (n + 3)4](ga5)»+3+". 

In these formulae, n represents the slope of the linear 
power spectrum at A:, i.e. n = d\n.P^(k) / d\n.k (with 
an additional spline interpolation as described in [14]), 
q = k/k^\ with /Cni defined by /c^iP/'(A:ni)/(27r^) = 1, and 
the function Qz{n) is given by 



4_2n 

1 -f- 2^+1 ' 



(21) 
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ai 


= 0.484, 


^2 


= 3.740, 


as = 


-0.849, 




= 0.392, 


as 


= 1.013, 


ae = 


-0.575, 


ar 


= 0.128, 


as 


= -0.722, 


ag = 


-0.926 



The best-fit values for the free parameters ai — ag were 
found by simulations [14] 



(22) 



III. PRIMORDIAL NON-GAUSSIANITY 

Additional to the contribution Bf^^^ from nonlinear 
gravity, the matter bispectrum can have primordial con- 
tributions B^^^ from inflation or some other early uni- 
verse model such as cosmic defects. While the simple 
model of single field slow roll inflation gives only a small 
primordial bispectrum, /nl ^ (9(10~^), other models can 
yield large non-Gaussianities with /nl > 1 (see e.g. [40- 
44] for reviews). Such models can be distinguished if 
they induce different bispectrum shapes, i.e. different de- 
pendencies of the bispectra on the momenta /ci , A:2 , ks 
as illustrated in Fig. 3. However, before reviewing pri- 
mordial bispectrum shapes we describe how primordial 
non-Gaussianity changes the dark matter bispectrum. 



Primordial contribution to the matter bispectrum 

Let us assume that an inflationary model produces a 
primordial bispectrum B^ with nonlinear amplitude /nl , 
i.e.2 

(^(ki)^(k2)^(k3)) = {27rfSD{^iki)hj^B^{k,,k2,ks). 

(23) 

From the linear Poisson equation (2) we see that the pri- 
mordial contribution to the matter bispectrum is given 
at leading order by 

5f"(A;i,fe,fe)U = 

M{ki,z)M{k2, z)M{ks, z)hi^B^{kiM. ks), (24) 

which is valid on large scales. 

A simple improvement to the tree level shape which 
incorporates some loop corrections can be obtained with 
the nonlinear power spectrum Ps: 



oprim/. u / P5iki)Ps{k2)P5{k^) p (u u u \ 

Bl^i^[ki,k2,k^) = 1/ ry n.\Ty n. \ TD n. ^ B^{ki, k2, k^) , 



P^{ki)P^{k2)P^{ks) 



(25) 



^ Non-linear corrections to the mapping between IZ and which 
are not taken into account here, induce an additional bispectrum 
of which is however smaller than the primordial contributions 
considered here [42]. 



where P^{ki) refers to the primordial power spectrum at 
some early time in matter domination. As will be demon- 
strated later in this paper, this shape can be used to 
obtain simple fitting formulae for the primordial contri- 
bution to the matter bispectrum. A more systematic but 
also more cumbersome approach is to include loop correc- 
tions, which become important on small scales [11, 12]. 
We shall calculate more of these contributions in a subse- 
quent paper and test their correspondence to the AT-body 
simulations [35]. 

We note that, while the time dependence of the tree 
level gravitational bispectrum (7) is given by Bf^^^ oc 
I^^(z), the tree level primordial contribution (24) only 
grows like B^™ oc D^{z), implying that it is easier to 
extract the primordial contribution to the dark matter 
bispectrum at early times. The simulations and fitting 
formulae discussed later will show that the gravitational 
bispectrum also grows faster than the primordial contri- 
bution in the strongly nonlinear regime (also by a factor 
of roughly D{z)). 



Non-Gaussianity as a halo model time-shift 

At sufficiently small scales the perturbative treatment 
breaks down and simulations and phenomenological mod- 
els must be used. The phenomenological halo model 
prediction for local non-Gaussian initial conditions was 
computed in [38] by incorporating modified expressions 
for the halo model ingredients in presence of local non- 
Gaussianity. While first tests by [38] demonstrate that 
this approach works well for some one- dimensional bis- 
pectrum slices in the mildly nonlinear regime, k < 
0.3/i/Mpc and z < 1, comparisons to simulations in the 
strongly nonlinear regime have not been undertaken to 
date. 

An alternative phenomenological model we propose 
here is to note that the primordial bispectrum can also 
contribute to the halo model bispectrum as an initial off- 
set or time-shift. As we have seen the 1-halo model is 
highly correlated with the simple 'constant' bispectrum 
we described earlier (14). The key point is that this con- 
stant halo contribution grows much more rapidly than 
both the tree-level primordial bispectrum (24) or the 
tree- level gravitational bispectrum (7). Consequently, if 
the primordial bispectrum has a significant positive (or 
negative) constant component, then this can act as an 
initial condition for the halo bispectrum; the faster halo 
amplification growth will start earlier (or later) by a time 
or redshift offset Az. The constant part of the primor- 
dial signal, then, will participate in the halo bispectrum 
growth and can be amplified much faster than expected 
from the tree-level result (24) or even with loop correc- 
tions (25). If this physical picture is correct, then the 
primordial constant contribution can be described per- 
turbatively around the constant halo bispectrum (14) by 



^ = 30 
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Figure 3. Illustration of separable expansion (44) for some theoretical tree level bispectra. Each black cell of 3 plots contains from left 
to right: Theoretical tree level prediction (7) or (24), separable expansion of the tree level theory (44) and estimated bispectrum from 
A/'-body simulations (56) (for simulations G512b, LoclO, EqlOO, OrthlOO and Flat 10 from top to bottom). The plots show signal to noise 
weighted bispectra (as on the left hand side of (44)) on the tetrapyd domain in Fig. 1, evaluated at redshift z = 30 for 512^ particles in 
a box with L = 1600Mpc//i. The opaqueness of the points reflects the absolute value of the signal to noise weighted bispectrum (values 
close to are completely transparent). Colors and transparency are consistent within each black cell but differ across different black cells. 
The plot axes are /ci, /c2, /cs < 0.5/i/Mpc. For better visibility we do not plot points with ki -\- k2 + ks > 2/cmax (corresponding to the green 
region in Fig. 1). 
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expanding the growth factor 



a time-shift for the more general 'constant' model (16) 



D'"'{z + Az) ^ D'"^{z) + nhD">^-\z)'^^^^Az. (26) 

dz 

We determine Az by matching the projection of the total 
measured bispectrum Bs on the constant shape ^o^nst 
in the Gaussian and non- Gaussian simulations: 



B 



( /const \ 
UNL jGa 



,{z + Az) 



(/nl'^^)ng(^)- 



(27) 



From simulations, the time shift Az and the correspond- 
ing shift of the growth function, AD = Az, are both 
found to vary over time at most by a factor of 2 for red- 
shifts 10 > 2: > 1 for local, equilateral and flattened ini- 
tial conditions (to the extent to which this can be tested 
by linearly interpolating the limited number of output 
redshifts at which the bispectrum was measured). Pro- 
vided AD is time-independent this perturbation of the 
simple 'constant' model (14) means we should be able to 
model the halo contribution from the primordial pertur- 
bation as 



BTZAkl,k2, ks) ^ C2 {Z) (fci + fc2 + , 



(28) 



where we expect n^^^ = Uh — I from (26) and again 
ly ^ —1.7. The fitting parameter C2 will be related to 
the correlation between the primordial shape^ and the 
constant model at the time at which halos form for the 
length scale under consideration. Different non-Gaussian 
bispectrum models should show consistent behaviour in 
the nonlinear regime depending on the relative magni- 
tude of their constant component. We shall define these 
quantities more precisely after introducing the bispec- 
trum shape correlator in the next section, but (28) will 
be an important component in our later fitting formulae 
for primordial non-Gaussianity. 

The simple power law in the linear growth function 
D{z) used to model the time dependence of the constant 
halo bispectrum can of course be extended to more gen- 
eral functions of time (16), whose time derivative would 
then enter in the expansion (26). The results presented 
later show that the overall normalisation of the simple 
fits can be improved by a more general modeling of the 
time dependence (reducing the overall growth rate in the 
strongly nonlinear regime). In this case, we can consider 



^ We refer here to the excess bispectrum compared to Gaussian 
initial conditions as measured in A/'-body simulations at times 
when the constant contribution to the bispectrum is not negligi- 
ble compared to the partially loop-corrected tree level contribu- 
tion (25). Therefore C2 is a fitting parameter to be determined 
from A/'-body simulations. 



5, const 

{ki,k2, ks)]^ = -Bf,const(^l' ^2, k3)\z+Az 

-BLZAkuk2,ks)U 



dT{k^ z, Zi) 
dz 



B^^Uk,zOAz 



(29) 



where we neglect higher order terms assuming them to 
be subdominant. We can determine the time-shift Az by 
determining the magnitude of the constant part of the 
primordial bispectrum at z = Zi. 

While more sophisticated models of the time evolu- 
tion are left for future work, we note that the correla- 
tion with the measured bispectrum shape cannot be im- 
proved much, because simulations are so well described 
by a combination of the 'constant' shape (28) and the 
modified tree-level gravitational shape (25) used in the 
fitting formulae presented later. 



Primordial bispectrum shapes from inflation 

Local shape 

The fiducial /nl model of primordial non-Gaussianity 
is the local model [45], which is described in this way be- 
cause it can be generated simply by squaring a Gaussian 
field in real space, 

$(x) = $g(x) + /nl[*^(x) - <$^>], (30) 

where {^q) denotes an average over x space and ensures 
that the average perturbation is zero. The resulting bis- 
pectrum takes the form 

B^i^{ki,k2, ks) = 2 [P^{ki)P^{k2) + 2 perms] , (31) 

which is illustrated in Fig. 3 and peaks at squeezed tri- 
angle configurations, where one wavenumber is much 
smaller than the other two. Multiple field infiation mod- 
els are one potential source of this shape (for a review 
see [40]). If a bispectrum signal in the squeezed limit is 
detected, then this will rule out all canonical single field 
models of inflation [46-48] . 



Equilateral shape 

Higher derivative operators in the inflationary action, 
arising e.g. in DBI inflation [49] and in effective fleld 
theory approaches [48, 50], produce a shape that can 
be approximated by the separable equilateral template 
[24, 40, 51] 

= 6[ - (P^{ki)P^{k2) + 2 perms) 

-2{P^{ki)P^{k2)P^{k^)f'^ (32) 

+ {Pl'\ki)Pl'\k2)P^{k^) + 5 perms)] , 
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which peaks in the equilateral limit, ki = k2 = ks. For 
equilateral triangles with ki = k2 = ks summing up the 
corresponding three plane waves, ^(kj)e*^J^, gives 
filamentary over densities, i.e. overdense cylinders along 
the direction perpendicular to the plane of the triangle 
of (ki,k2,k3), surrounded by underdensities (as moti- 
vated in Ref. [52]). A primordial equilateral shape from 
higher derivatives, therefore, must be distinguished care- 
fully from the stronger equilateral and (highly correlated) 
constant contribution produced at late times by nonlin- 
ear gravitational collapse and the emergence of filamen- 
tary and point-like structures (see e.g. [31, 33]). We will 
confirm this expectation by measuring the bispectrum in 
A^-body simulations. 

Orthogonal shape 

Another shape that can arise from single field inflation 
is the orthogonal shape [53] , which is roughly orthogonal 
to the equilateral and the local shape and peaks (with op- 
posite sign) for both equilateral and flattened/elongated 
triangle configurations {ks = ki-\- k2). It can be approx- 
imated by the separable template [53] 

= 6[3{Pl^^{ki)Pl^\k2)P^{k3) + 5 perms) 

- Is'i' -8{P^{k,)P^{k2)P^{ks)f^']. 

(33) 

Being orthogonal to the equilateral shape, the constant 
mode is relatively suppressed, a fact which will be im- 
portant in later discussions. We note that all the shapes 
above (31)-(33) are written explicitly as the sum of sep- 
arable functions of the wavenumbers /ci , A:2 , k^. 

Flattened shape 

A non-Bunch-Davies vacuum leads to shapes which 
peak for flattened/elongated {ki -\- k2 = ks) and folded 
{ki = 2/c2 = 2ks) triangle configurations'^ [40, 54- 
56]. This flattened shape depends on terms such as 
1 / {ki -\- k2 — ks) [54]. Therefore it provides an impor- 
tant example of a bispectrum which is inherently non- 
separable, so it is computationally difficult to generate 
initial conditions for this shape and to extract its ampli- 
tude from data. We will use an expansion in separable 
modes to solve both problems [7, 8, 26, 27]. Of course, 
there is a divergence for elongated triangles, ks = /ci +A:2, 
which must be removed with a physically motivated cut- 
off. For definiteness we use the same shape and cutoff 



Sometimes other names are used for these triangles, we follow 
[34]. 



as was used for the CMB in [2, 7], setting B = for 
ki-\-k2 — ks < 0.03(/ci -\-k2-\-ks) and permutations, and 
then smoothing on ki -\- k2 -\- ks = const, slices with a 
Gaussian filter with FWHM of 0.03/(A:i + + ks). 

It is worth noting that the equilateral and orthogo- 
nal templates are separable approximations of inflation- 
ary shapes, with good agreement where the primordial 
bispectrum signal peaks. However late time observables 
may be sensitive to suppressed triangle configurations 
for which template and physical shape may differ signifi- 
cantly (e.g. scale dependent halo bias mainly depends on 
the squeezed limit of the bispectrum [24, 57]). The meth- 
ods described below enable us to simulate and estimate 
physical non-separable bispectra without using such sep- 
arable templates. However we do not expect significant 
differences for the dark matter bispectrum, because its 
peaks are well approximated by the separable templates. 
For this reason and for easier comparability with other 
simulations, we use the equilateral and orthogonal tem- 
plates instead of the physical shapes. Future work, es- 
pecially on the power spectrum and bispectrum of halos, 
will instead focus on the physical shapes. 



IV. BISPECTRUM ESTIMATION 
METHODOLOGY 

/nl estimator 

If our theoretical model is that the density perturba- 
tion S has power spectrum Ps{k) and bispectrum f^Bf^, 
then the maximum likelihood estimator for the amplitude 
of this bispectrum in the limit of weak non-Gaussianity 
is given by [26, 27]^ 

Bf'i.ki, fc2, k3)[dk,5kji,^ - 3(Jki^k2)43] /o4^ 
Ps{ki)Ps{k2)Ps{h) ^ ' 

where 5 is the observed density perturbation. If this has 
a bispectrum Bf^^^ i.e. 

(^k.^k.^ks) = (27r)^^z^ (E,k,) Bf%kuk2, ks), (35) 



^ In the denominator we have assumed that (SS) is diagonal which 
is valid for a statistically homogeneous density perturbation. The 
linear term {SS)S in the numerator is written out for completeness 
but it is not used in our simulations because it vanishes for a 
statistically homogeneous field (because S\^=q = 0). 
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then the expectation value of (34) is given by^ 



1 V 

Nth TT jvb 



dVk kik2k3 



Psiki)Ps{k2)Ps{k3) 



(36) 



where dVk = dkidk2dks, Vb is the tetrahedral domain 
aUowed by the triangle condition on the wavenumbers 
ki^ and ^ is a volume factor given by ^ = {27r)^6D{0) = 
L^. Demanding (f^^) 
normalisation such that 



1 for 5f 



Bf^ fixes the 



th 



^ Jvb 



dVk 



klk2k3[Bf{ki,k2,k3)]' 

Ps{ki)Ps{k2)Ps{k3) 



(37) 



Shape and size comparisons 

Eq. (36) motivates the definition of a scalar product 
between two bispectrum shapes [26, 58], 



{Bi,B,) = - [ dVu 



kik2k3Bi{ki,k2j k^)Bj{ki,k2, ks) 
P5{h)Ps{k2)P5{ks) 



(38) 

which can be normalised to a number between —1 and 1 
by defining the cosine or shape correlation 



C{Bi,Bj) 



{Bi , Bj 



y^{Bi,Bi) (Bj.Bj) 



(39) 



If two theoretical bispectra Bi and B2 have a small shape 
correlation, \C{Bi, B2)\ <C 1, the optimal estimator for 
/nl ^^^^ perform badly in detecting and vice versa. 
Therefore the shape correlation is a useful measure of the 
similarity of bispectrum shapes. 

To measure the total integrated size of a bispectrum 
we define the squared norm as the cumulative signal to 
noise squared of the bispectrum [2, 51], 



m 



V 



(27r)a 



SB B'^{ki,k2,kz) 
6 var 5(^1,^2,^3) 



{B,B) 
6(27r)3' 
(40) 



which equals the quantity {S/N)\ in [17]. The symmetry 
factor 5b is 6 if /ci = = , 2 if only two ki equal each 
other and 1 if all ki are different from each other. To 
obtain the expression on the right hand side of (40) we 
assumed that different Fourier modes are uncorrelated, 



^ Performing the angular integrals shows for arbitrary F{ki , /c2 , /cs) 
which corrects for a factor of (27r)^ missing in [26, 27]. 



(^k^k') 0^ ^^(k + k'), and we only took the Gaussian 
contribution to the bispectrum noise into account, i.e. 



Ydjc B{ki,k2,k^) = (27r) 



sSBP6{h)Ps{k2)P6{ks) 



87r^ kik2ks 



(41) 



It is also convenient to normalise the total integrated 
bispectrum size with respect to the linearly evolved local 
bispectrum defined in (31) and (24) (see [2] for a similar 
quantity in the CMB context): 
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NL 



prim,loc i 



\B 



(42) 



We use the linearly evolved local bispectrum without 
loop corrections to obtain an expression which can be 
easily evaluated. 

Using definitions (39) and (40), the expectation value 
(36) can be rewritten in the intuitive form 



(/nl) 



C{Bf,Bf^) 



\\Bn 



(43) 



Therefore the bispectrum amplitude /nl can be inter- 
preted as the projection of the observed bispectrum shape 
on the theoretical bispectrum shape, which is given by 
the cosine of the shapes times the ratio of their norms. 



Separable mode expansion for fast /nl estimation 

Evaluating the estimator (34) is computationally very 
expensive because it requires 0{N^) operations for a grid 
with TV points per dimension (typically N = O(IO^)). 
For an efficient evaluation we expand the signal to noise 
weighted theoretical bispectrum in the separable form 



VhhhBf'{ki,k2,ks) 
^yPs{ki)Ps{k2)Ps{k3) 



^ a^q{r{ki)qs{k2)qt}{k3) , 
(44) 



n=0 



where {rst} denotes a symmetrisation over the indices, 
and a partial ordering has been introduced to enumerate 
the modes (for a fuller description see [8]). The expres- 
sions Qr could be any set of independent one-dimensional 
basis functions. For definiteness we choose Qr to be poly- 
nomials of order r defined in [8] and choose the slice or- 
dering defined in Eq. (58) of [8]. We truncate the modal 
expansion (44) after rimax = O{50) modes, which pro- 
vides an accurate representation for the shapes we are 
considering as illustrated in Fig. 3, which compares the- 
oretical bispectra on the left with truncated expansions 



^ We need an additional factor of (27r)^ compared to [17] because 

^them = (27r)3Pu3 and Bthem = (27r)35us. 
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in the centre (and bispectra measured in N-body simula- 
tions on the right). More quantitative discussions of the 
error induced by the truncation of the expansions wiU be 
provided later in terms of shape correlations. 

The separable expansion (44) allows us to write the 
estimator and its expectation value in the form 



/i 



■th 
NL 



{2nl 
Nth 



-3(M{,(x)M,(x))M,}(x)], (45) 
ifNl.) —T.'^^^^lnm, (46) 



where 



V 



/ dVkQn{kl,k2,h)Qm{kl,k2,k3), 



,^ / ^ f ^ • 



IVb 



f ^kPM 



(47) 



and 



Qn{k, k\ k") = qir{k)qs{k')qt}{k") . (48) 



This form of the estimator can be evaluated efficiently 
because it involves only Fourier transforms and one 
three-dimensional integral over position space, leading to 
0{N^) operations. 

To motivate the choice of the prefactor in the 
expansion (44) let us transform to basis functions 
Rn{ki,k2,ks) that are orthonormal on the tetrapyd do- 
main, (i^n, ^m)noweight = <^nm, with rcspcct to the Un- 
weighted scalar product 



V 



J I noweight 



Vb 



dVk Bi{ki,k2, ks)Bj{ki,k2, ks). 

(49) 



From J2n o^nQn = Z)„ ct^-Rn and Eq. (44) we find that 
Bfik^MM) = J2^nB^{kuk2,k3), (50) 

n 

where the contributions 



B^{kuk2M 



I Ps{ki)Ps{k2)Ps{k3) 



k\k2kz 



Rn{ki,k2,kz) 



(51) 



to the full bispectrum are orthonormal with respect to 
(38). Thus the expansion coefficients measure the size 
of contributions to the bispectrum which are orthonor- 
mal with respect to the signal to noise weighted scalar 
product (38). 

As a side remark, note that in contrast to the signal 
to noise weighted scalar product (38) induced by the es- 
timator expectation value, the weight of the scalar prod- 
uct (49) is time-independent and does not depend on the 



range of scales under consideration. Therefore we can use 
the same set of orthonormal polynomials Rn at any time 
and for all length scales. If the expansion (44) does not 
converge as well as for the applications considered here, 
different (separable) scalar product weights can be used 
to orthonormalise the RnS. 



Fast modal bispectrum estimator 

Importantly the separable mode expansion allows us 
not only to measure efficiently for any given the- 
oretical bispectrum, but it also allows us to recon- 
struct the full bispectrum in a model-independent man- 
ner by measuring a set of independent /nl's and sum- 
ming up the corresponding contributions to the bispec- 
trum [26, 27]. To see this, note that the transformation 
= J2rn ^nmQm implies = J2p KnOL^ and from 

(i?n,^m) noweight = ^nm WC get 7 = A"^(A"^)^. Plug- 
ging this into (45) gives 



/nl 



1 



R 



(52) 



with 



E(27r)3^A„„ j S^[MrMsMt-Z{M{rMs)Mt}\, 

n ■' 

(53) 



where the index n labels the combination r, s, t. If Bf^ 
{Bf^) then (46) becomes 



Equations (52) and (54) imply 

(/?^) = « 



(54) 



(55) 



Therefore we can estimate the full bispectrum with [26, 
27] 



Bs{kiMM)= £ l3nB^{ki,k2,k3), 



n=0 



(56) 



where B^ is given by (51). Intuitively the coefficients 
measure the components of the orthonormal basis 
bispectra B^ in the data and Bs measures the projection 
of the bispectrum on the subspace of bispectra spanned 
by the B^. 

This approach is extremely efficient with the calcula- 
tion of each coefficient being equivalent to a single 
3D integration over products of fast Fourier transforms. 
While theoretically a complete basis would require 
modes in practice far fewer modes (0(50)) are necessary 
to reconstruct the bispectrum accurately. 
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The primordial contribution to the matter bispectrum 
win be extracted by measuring the difference of the bis- 
pectrum between non-Gaussian and Gaussian initial con- 
ditions, 



^NG = ^2 ~ (/^n) Gauss] 



(57) 



n=0 



As will be discussed in more detail later, the difference is 
computed seed by seed to reduce error bars. 

The time-dependence of the expansion coefficients 
and the (3^ coefficients can be read off from (44) to be 
(<^n)grav D{z) foT the gravitational bispectrum and 
(q^^)ng oc D(z)^ for primordial bispectra, assuming tree 
level perturbation theory where Ps oc D'^{z). In practice 
we must use the nonlinear power in (44) and (47) im- 
plying slightly different time dependences for in the 
nonlinear regime (see e.g. the middle panel of Fig. 14a, 
which will be discussed in more detail later). 



Fast modal bispectrum correlations 

To analyse the reconstructed bispectrum and to control 
the accuracy of the separable method, the following shape 
correlations are computed using (39): 




\ n 



th Dth\ 



(58) 



A^m\ mJ A^pXr'p / sim 
C/3,th =C [{Bs)^ira,Bf'^ = Ci3^aCa,th- (60) 

Here (Bs) sim. denotes the average of (56) over indepen- 
dent simulations and we truncate all sums appearing 
in the above expressions and in (56) after rimax modes. 
The first shape correlation measures how well the sep- 
arable expansion (44) of the theoretical bispectrum ap- 
proximates the theoretical bispectrum. The second shape 
correlation quantifies how well the estimated bispectrum 
{Bs) sim agrees with the separable expansion of the theo- 
retical bispectrum. The product of these two shape cor- 
relations gives the shape correlation between the recon- 
structed bispectrum and the theoretical bispectrum. 

Cumulative measures of non-Gaussianity 



= , (59) 



noise of the average reconstructed bispectrum, which can 
be expressed in terms of the measured average {(3^) sim 
as^ 



{Bs)s 



R\2 
n I sim 



6(27r)3 



ZTisim 



\Enm 



R\2 
sim 



-A^^prim,loc 



(61) 



The norm was defined in (40) using the cumulative sig- 
nal to noise squared of the bispectrum. As a consistency 
check we check if (61) also measures the cumulative signal 
to noise squared of the f3^ coefficients. The theoretical 
covariance of involves the 6-point function of the den- 
sity perturbation. Taking only the Gaussian contribution 
into account and assuming that different Fourier modes 
are uncorrelated we find 



(62) 



This predicts ct^r = Y^6(27r)^ for Gaussian simulations 
and therefore confirms that (61) measures the cumulative 
signal to noise of the reconstructed bispectrum. We con- 
firmed the prediction for a^R quantitatively in Gaussian 
simulations, but we have not assessed non- Gaussian cor- 
rections to the predicted noise (see e.g. [59]). Note that 
our plots show la sample standard deviations obtained 
by running realisations with different random number 
seeds. 



Towards experimental setups 

For the bispectrum measurements presented below we 
will drop the second term in the square brackets in (53) , 
i.e. we assume that different modes of the density pertur- 
bation are uncorrelated, (^k^k') (^^(k -h k'), which is 
true for a statistically homogeneous field. In experimen- 
tal setups, which are not considered here, off-diagonal 
mode couplings from inhomogeneous noise can be incor- 
porated in (53) with the linear term (MM)M, without 
affecting the efficiency of the bispectrum estimator (56). 
Ref. [59] explores how modifications to the denominator 
of the estimator (34) in presence of off-diagonal mode 
couplings can be incorporated efficiently with modal ex- 
pansions and an implementation for CMB experiments 
has been developed successfully. 



Comparison to other bispectrum estimators 

An alternative way to reconstruct the bispectrum is 
to estimate the three-point correlation function directly. 



Additionally to these shape correlations we will also 
compare the projections (52) of the measured bispectrum 
on theory bispectra, as well as the cumulative signal to 



Alternatively using (X^^(/3^)^)sim would add undesirable con- 
tributions from the variances of the /3^, which become relevant 
if the measurements are not signal-dominated. 
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taking a subset of all possible triangle configurations that 
yields the desired accuracy of the bispectrum [12, 14, 25]. 
In contrast our estimator takes all possible triangle con- 
figurations into account in a very efficient way, so that 
the bispectrum can be measured on the complete range 
of scales that were included in the A/'-body simulation. 
For example, with our method the extraction of the bis- 
pectrum of a 1024^ grid takes about one hour on 6 cores, 
which is only a small fraction of the time required to 
run an A/'-body simulation of this size. Moreover our es- 
timator compresses the information about the shape of 
the bispectrum, which is a function of the whole three- 
dimensional tetrahedral domain allowed by the triangle 
condition, to nmax numbers where nmax is the num- 
ber of basis functions used in the expansion (44). This 
simplifies further processing of the bispectrum, e.g. for 
comparisons with theoretical models and for obtaining 
fitting formulae. Once is measured from the data, we 
can not only obtain the full bispectrum with (56), but we 
can also calculate full three-dimensional shape correla- 
tions, the cumulative signal to noise and the nonlinearity 
parameter associated to any theoretical bispectrum 
Bf^ without additional computational cost, using (59), 
(61) and (52), respectively. In presence of inhomogeneous 
noise we can in principle include off-diagonal covariance 
elements {S\^S\^f ) in a straightforward way. Finally our ap- 
proach allows us to estimate the trispectrum efficiently 
[26], which has been shown for a class of trispectrum 
shapes in [27] (see also [4, 9]), but we leave the applica- 
tion to A^-body simulations for future work. 

Bispectrum visualisation 
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Figure 4. Bispectrum weights [Ps (ki) Ps (k2) + 2 perms] ^ 
(top) and kik2k3 / [Ps (ki) Ps (k2) Ps (ks)] (bottom) evaluated with 
CAMB [28] at redshift z = 30 on shces with ki+k2 + k3 = 1/i/Mpc. 



Often in the literature the bispectrum is visualised 
by plotting one- or two-dimensional slices through the 
tetrapyd shown in Fig. 1. E.g. the plots in [58] corre- 
spond to two-dimensional slices through the tetrapyd ob- 
tained by varying k2 and ks at fixed ki and removing the 
ks > k2 part of the slice by symmetry. Some plots in [7] 
and in this paper show two-dimensional tetrapyd slices 
with ki -\- k2 -\- ks = const. While such slices are sufficient 
for scale- invariant primordial bispectra, late time bispec- 
tra typically have different triangle dependences at dif- 
ferent overall scales ki-\-k2-\-ks. This motivates plotting 
late time bispectra on the full three-dimensional tetrapyd 
instead of plotting particular slices through the tetrapyd. 

Instead of the unweighted bispectrum the so-called re- 
duced bispectrum Q = B§ / {P§{ki)Ps{k2) + 2perms) is 
often shown, because it reduces the dynamical range of 
the bispectrum and for the tree level gravitational contri- 
bution it is independent of time, overall scale and power 
spectrum normalisation and almost independent of cos- 
mology [31]. Instead we will plot the signal to noise 
weighted bispectrum kik2ks / {P5{ki)P5{k2)P5{ks))Bs 
because then the product of the shown functions gives 



the scalar product defined in (38) [58], visual similarity 
indicates a high shape correlation (39) and the dynam- 
ical range of the unweighted bispectrum is also greatly 
reduced, which is of advantage for plotting purposes. For 
linearly evolved primordial bispectra (24), the signal to 
noise weighted bispectrum is time-independent in the lin- 
ear regime. 

Qualitatively the two weights are quite similar in the 
regime relevant for most plots of this paper, see Fig. 4. 
However they differ in the squeezed limit, because for 
decreasing /ci, Ps{ki)Ps{k2) turns over at ki = /Ceq^ 
whereas y^P6{ki)Ps{k2)P6{k3)/{k 

1^2 ^3) turns constant 
where P{ki)/ki turns constant, which is on somewhat 
larger scales than /ceq. While it is not straightforward to 
deduce the squeezed limit of an unweighted bispectrum 
from a plot of its weighted form, one can compare plots of 
different bispectra if they are weighted in the same way 
to deduce the relative behavior in the squeezed limit. 



V. SIMULATION SETUP, INITIAL CONDITIONS 
AND VALIDATION 

A/^-body simulations setup 

We use the separable mode expansion method de- 
scribed in [26, 27] to generate reahsations of the initial 
primordial potential ^ = —37^/5 during matter domi- 
nation with the desired primordial power spectrum and 
bispectrum. For 512^ particles this takes about 10 min- 
utes per seed on one core and works for separable as well 
as non-separable bispectra, therefore being more general 
than other proposed methods [23, 25]. From <l> we calcu- 
late the linear density perturbation S at the initial red- 
shift of the simulation with the Poisson equation (2). We 
then use this initial density perturbation to displace the 
initial particles from an unperturbed distribution with 
the 2LPT method [60, 61], which also determines the ini- 
tial particle velocities. Then the A^-body code Gadget-3 
[62, 63] simulates the time evolution until today and we 
use a cloud in cell scheme to calculate the density pertur- 
bation 6 of the particle distribution on a grid at different 
redshifts. After deconvolving 6 with the cloud in cell 
kernel we compute the power spectrum Ps and the coef- 
ficients f3^ from (53) using nmax = 50 modes. Finally we 
reconstruct the full bispectrum with (56) and calculate 

its norm Fnl (61) as well as its shape correlation C^^a 
(59) with theoretical bispectra and its nonlinear ampli- 
tude (52). 

Table I lists the parameters of the A/'-body simula- 
tions that were performed in this work. All simula- 
tions assume a flat ACDM model with the WMAP- 
7 [1] parameters n^h'^ = 0.0226, r^c/^^ = 0.11, f^A = 
0.734, /i = 0.71, r = 0.088, A2^(A:o) = 2.43 x 10"^ and 
ns{ko) = 0.963, where ko = 0.002Mpc"^ Note that 
the primordial power spectrum is given by P^{k) = 
(9/25) (27rVfc3) A2,(fco)(fc/fco)"'-i. 



Regular grid vs glass initial conditions 

For the unperturbed particle distribution, from which 
initial particles are displaced using the 2LPT method, we 
either use a regular grid or a glass configuration, which 
is obtained by placing particles randomly in the box and 
then evolving them with the sign of gravity flipped in 
Gadget-3 [63, 64]. A disadvantage of glass initial con- 
ditions is that they require generating a glass and spec- 
ifying when the glass configuration has converged (see 
e.g. [65] for a discussion). In contrast, regular grid ini- 
tial conditions are easier to set up but break statistical 
isotropy, which may have implications for structure for- 
mation. We will compare the diff"erent initial condition 
setups first for Gaussian and then non- Gaussian initial 
conditions. 
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Name 


NG 


/nl 




Np 






Nr 


glass 




shape 














G512g 


- 


- 


1600 


512 


49 


156 


3 


yes 


G512 


- 


- 


1600 


512 


49 


156 


3 


no 




- 


- 


{400, 100} 


512 


49 


{39,9.8} 


3 


no 


G768 


- 


- 


2400 


768 


19 


90 


3 


no 


G1024 


- 


- 


1875 


1024 


19 


40 


2 


no 


LoclOg 


local 


10 


1600 


512 


49 


156 


3 


yes 


LoclO 


local 


10 


1600 


512 


49 


156 


3 


no 


LoclOi^^ 


local 


10 


{400, 100} 


512 


49 


{39,9.8} 


3 


no 


LoclO" 


local 


-10 


1600 


512 


49 


156 


3 


no 


Loc20 


local 


20 


1600 


512 


49 


156 


3 


no 


Loc50 


local 


50 


1600 


512 


49 


156 


3 


no 


EqlOOg 


equil 


100 


1600 


512 


49 


156 


3 


yes 


EqlOO 


equil 


100 


1600 


512 


49 


156 


3 


no 


EqlOOi^^ 


equil 


100 


{400, 100} 


512 


49 


{39,9.8} 


3 


no 


EqlOO" 


equil 


-100 


1600 


512 


49 


156 


3 


no 


OrthlOOg 


orth 


100 


1600 


512 


49 


156 


3 


yes 


OrthlOO 


orth 


100 


1600 


512 


49 


156 


3 


no 


OrthlOOtJg 


orth 


100 


400 


512 


49 


39 


3 


no 


OrthlOO" 


orth 


-100 


1600 


512 


49 


156 


3 


no 


FlatlO 


flat 


10 


1600 


512 


49 


156 


3 


no 


FlatlOtJg 


flat 


100 


400 


512 


49 


39 


3 


no 



Table 1. Parameters of A/'-body simulations: Non-linearity pa- 
rameter /nl, box size L, number of particles per dimension 
Np, initial redshift of the simulations Zi, softening length Lg 
and number of realisations (i.e. random seeds) Nr for each 
parameter set. 'glass' indicates if the initial particles were 
displaced from a regular grid or from a glass conflguration. 
Initial conditions for non-local non-Gaussian simulations were 
generated with the separable method described in [26, 27]. 
All simulations use 2LPT [60, 61] to get the initial particle 
distribution, which is then evolved with Gadget-3 [62, 63]. 



Measured bispectra for Gaussian A/'-body simulations 
will be shown in Fig. 14 in Section VI. Simulations 
started from glass and regular grid initial conditions are 
shown in Fig. 14a and Fig. 14b, respectively. At early 
times regular grid initial conditions produce a bispec- 
trum which is more than three times as large as the tree 
level prediction. This signiflcant spurious bispectrum is 
not present for glass initial conditions and must therefore 
be caused by the breaking of isotropy induced by the reg- 
ular grid. Fig. 14b shows that the spurious bispectrum 
decays with time as the regular grid structure is washed 
out by the growing gravitational perturbations. At late 
times, z < 3^ the bispectra from regular grid and glass 
initial conditions differ by at most 10% in their cumula- 
tive signal to noise and by less than 0.5% in their shape 
correlation to the tree level prediction, see Fig. 14c. 

We conclude that at the 10% accuracy level (for Gaus- 
sian simulations with L = 1600Mpc//i, 512^ particles and 
^max = 0.5/i/Mpc), both glass and regular grid initial 
conditions can be used to extract the gravitational bis- 
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Figure 5. Impact of glass initial conditions on measured non-Gaussian bispectra in simulations LoclO, EqlOO and OrthlOO. We plot the 
ratio to simulations with regular grid initial conditions. The top panel contains error bars obtained by calculating the sample standard 
deviation of the ratio for each seed. We do not show error bars in the other panels because they depend on J^n (/^i?)sim therefore 
more difficult to estimate. However they should be similar to the errors in the top panel because of (43). The large deviations in the 
bottom panel are mainly due to the different gravitational bispectra for glass and regular grid initial conditions as shown in Fig. 14c. 



pectrum from Gaussian simulations as long as the bispec- 
trum is measured at low redshifts when the regular grid 
is washed out. The difference between regular grid and 
glass initial conditions decreases when reducing the cut- 
off /Cmax used for the bispectrum estimation because then 
structures of the size of the grid separation are smoothed 
out. 

For non-Gaussian simulations, Fig. 5 shows how glass 
initial conditions change the non-Gaussian excess bispec- 
trum B^Q (57) compared to regular grid initial condi- 
tions. At z < 10 the shape of ^ng is consistent between 
glass and regular grid initial conditions at the 1.5% level 
for the simulations LoclO, EqlOO and OrthlOO. However, 
the total integrated bispectrum Fnl differs by up to 7% 
for z < 10. At earlier times the deviations are some- 
what larger. Compared to the full gravitational bispec- 
trum for Gaussian initial conditions shown in Fig. 14a, 
the impact of glass initial conditions on ^ng is quite 
small, possibly because the spurious bispectrum due to 
the regular grid is partly subtracted out when calculating 
^NG = B — Bqs^^^^. 

Whether regular grid or glass initial conditions are rep- 



resenting the statistics of the non-Gaussian density per- 
turbation more successfully is not unambiguous. Regular 
grid initial conditions introduce a spurious contribution 
to the full bispectrum B at early times, which can be 
avoided using glass initial conditions. However glass ini- 
tial conditions represent the primordial contribution B^q 
slightly less accurately than regular grid initial conditions 
at the starting redshift of the simulation. However these 
effects impact B^q only at the (9(5%) level at z < 10, 
so that both glass or regular grid initial conditions can 
be used at this level of precision. As for Gaussian sim- 
ulations the two initial condition methods agree better 
when /Cmax is reduced (see Fig. 5) and grid discretization 
effects are minimized. 



Validation and convergence tests 

First we test the setup of the initial conditions and the 
A/'-body simulations by comparing the measured matter 
power spectrum with the power spectrum predicted by 
linear theory and by CAMB [28, 30] in Fig. 6. Next we 
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Figure 6. Measured matter power spectra in simulations Gf Jo 
(red), G4J0 (green), G512 (blue) and G512g (cyan) in compar- 
ison with power from linear perturbation theory (grey dashed) 
and CAMB (grey solid, [28, 30]). See Table I for param- 
eters. Curves from bottom to top correspond to redshifts 
z = 49,20,6,2,0. 



8000 
6000 
4000 
2000 





-6^ 

H n 



10 



20 30 
n 



40 



50 



Figure 7. Shot noise validation test for the separable bispec- 
trum estimator: Measured bispectrum coefficients (53) 
(blue, averaged over 3 seeds) for 512^ particles placed ran- 
domly in a L = 1600Mpc//i box, compared with the expan- 
sion coefficients (44) (green) for the expected pure shot 
noise bispectrum. The shape correlation is Cp^a = 0.9998. 
For better visibility the region under the curve is colored 
green. 



perform a simple test of the bispectrum estimator by dis- 
tributing particles randomly in a box and comparing the 
measured bispectrum with the pure shot noise bispec- 
trum 5f°* = L^/N^ [66], see Fig. 7. These two tests 
show that the basic setup of the simulations is correct 
and that the separable bispectrum estimator works well. 

To check convergence of the A^-body simulations and 
the bispectrum measurements we compare results for box 
sizes L = {1600,400, 100}Mpc//i with = 512^ parti- 
cles. For each simulation we measure the bispectrum 

= {1,0.5,0.25} X (^^^ 



imposing a sharp k filter on the density perturbation 5]^. 

The measured late time bispectra are shown in Fig. 8 
as functions of /cmax- Data points with the same /cmax 
measured in a low and high resolution simulation can 
differ in the plots because, for these data, /Cmin differs 
by a factor of 4 and, therefore, the number of modes 
for the bispectrum estimation differs by a factor of 4^ = 
64, which affects particularly the cumulative signal to 
noise ||^||. In contrast the quantity F^^ is normalised 
with respect to the linearly evolved local shape over the 
given range of scales and therefore takes differences in the 
number of modes into account. Thus the agreement of 
F^L (and C/3^q,) for the different simulations seen in Fig. 8 
shows that the A/'-body simulations have converged and 
the bispectrum estimator gives consistent results. 

It should be noted that the shown cumulative signal 
to noise takes only Gaussian noise into account, c.f. (41). 
This simplification is not expected to be valid in the 
strongly nonlinear regime so that the correct signal to 
noise may differ significantly. However to assess this is- 
sue in a robust manner we need to run more realisations 
of the A/'-body simulations or use larger boxes with more 
particles. 
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Error bars 



To obtain error bars for quantities based on the primor- 
dial contribution to the bispectrum, we use seed by seed 
subtracted coefficients — (/3^) Gauss )sim and calculate 
their sample standard deviations, which are then used for 
standard error propagation assuming uncorrelated j3^. 
The error bars therefore measure how much the primor- 
dial contribution to the matter bispectrum scatters for 
different seeds of the initial Gaussian field in the sim- 
ulations. Compared to calculating the difference between 
the average bispectra, (/3^)sim - ( (/^^ ) Gauss )sim, the seed 
by seed subtraction reduces the error bars, because the 
late time bispectra for Gaussian and non-Gaussian initial 
conditions are very correlated due to the presence of the 
large gravitational bispectrum in both cases. 

In real observations we do not know the realisation of 
our universe with perfectly Gaussian initial conditions 
and therefore error bars will be larger if we measure 
/^m - ((/^m) Gauss) Sim. We do uot discuss the interesting is- 
sue of observability of primordial non-Gaussianity here, 
because instead of dark matter bispectra this requires 
halo bispectra, which we leave for future work (but see 
[17] for the local shape). 
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Convergence tests: Dependence of measured bispectra on /cmax at z = for simulations with 512"^ particles and box sizes 



{1600,400, 100}Mpc//i (G512, G^ig^ ^^q^. 



LoclO, LoclO|4oo,ioo} 5 EqlOO and Eql00j-4QQ ^oo})- (^) '^^P P^^e/; Projection of measured 
bispectrum on the tree level prediction (7) (squares) and on the fitting formula from [14] (circles). Middle panels: Cumulative signal to 
noise of measured bispectrum (squares/circles), tree level prediction (7) (dashed lines) and fitting formula from [14] (solid lines). Bottom 
panel: Shape correlation of measured bispectrum with tree level prediction (squares) and fitting formula from [14] (circles), (b-c) Measured 
non-Gaussian bispectra ^ng (57) in simulations with local and equilateral initial conditions compared to the linearly evolved primordial 
bispectrum (24) (dashed lines). Note comments on the cumulative signal to noise on small scales in the main text. 



VI. GRAVITATIONAL BISPECTRUM RESULTS 
Gravitational collapse and bispectrum evolution 

The cosmic web simulated by A/'-body simulations has 
a complex filamentary structure, which is illustrated in 
Fig. 9a. The bispectrum of this dark matter distribu- 
tion is shown in Fig. 9b. In the following we will discuss 
how the bispectrum can be used as a quantitative tool to 
characterise the structures formed by gravitational col- 
lapse and to study modifications due to primordial non- 
Gaussianity. We will demonstrate the unbiasedness of 
our bispectrum estimator by comparison with perturba- 
tion theory for large scales at early times. 

We can develop a qualitative and visual understanding 
of the bispectrum by comparing it to the form of large- 
scale structures. Fig. 10 shows snapshots of the dark 
matter distribution at redshifts 2: = 4, 2 and on the 
left and the corresponding measured bispectrum signal 
on the right. It is apparent that the shape of the bispec- 



trum characterises the three-dimensional structures that 
have formed. The diffuse blob- and pancake-like struc- 
tures at early times, 2: = 4, correspond to a fiattened bis- 
pectrum, which peaks at the edges of the tetrapyd and is 
predicted by leading order perturbation theory. At later 
times we find that filaments and clusters induce a bis- 
pectrum with a relatively enhanced signal for equilateral 
triangle configurations, as expected from sections II and 
III. As illustrated in Fig. 11, the bispectrum signal shows 
a similar transition for varying overall scale ki-\-k2-\-ks at 
fixed time z = 0, which indicates self-similar behaviour. 



Comparison with leading order PT 

On large scales and at early times we expect the grav- 
itational bispectrum to agree with leading order pertur- 
bation theory (7). To test this Fig. 12a shows the f3^ 
coefficients (53) measured in the Gaussian A/'-body sim- 
ulations G512 with L = WOOMpc/h in comparison with 
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(a) (b) 

Figure 9. (a) Dark matter distribution in one of the G^qq simulations of 512^ particles in a box with L = 400Mpc//i at redshift z = 0. 
(b) Measured (signal to noise weighted) bispectrum in the range 0.016/i/Mpc < k < 2/i/Mpc, averaged over the simulation on the left and 
two additional seeds. The opaqueness of the plotted points reflects the absolute value of the bispectrum. 



the expansion coefficients (44) of the tree level gravity 
bispectrum (7). The corresponding theoretical and esti- 
mated bispectra (56) are plotted over the full tetrapyd in 
Fig. 3 for 2; = 30 and Fig. 13 for 2; = 2 and z = 0. At early 
times, z = 30, the bispectrum agrees with tree level per- 
turbation theory (7), which predicts that the signal to 
noise weighted bispectrum is large for folded and elon- 
gated configurations but is suppressed in the squeezed 
limit. However at 2; = for /c > 0.25/i/Mpc the TV-body 
bispectrum has an enhanced amplitude for elongated con- 
figurations and an additional equilateral contribution not 
captured by tree level perturbation theory, which breaks 
down on these scales as expected. 

To analyse the measured bispectrum more quantita- 
tively and to illustrate its time and scale dependence we 
show in Fig. 14 the measured amplitude /^i^^ of the grav- 
itational bispectrum, its cumulative signal to noise \\B\\ 
(40) and the cumulative signal to noise normalised to the 
local shape Fnl (42), as well as the shape correlation C 
(39) of the measured bispectra with the theoretical ex- 
pectation. The meaning of these quantities is illustrated 
in Fig. 14f. We plot them as functions of redshift z and 
use different colors for different k ranges used for the bis- 
pectrum measurements. Note that is the amplitude 
of the gravitational bispectrum (7), i.e. tree level pertur- 
bation theory predicts = 1 in this convention. 

Fig. 14a and Fig. 14b show that for 0.0039/i/Mpc < 
k < 0.5/i/Mpc at z = the tree level bispectrum under- 
predicts the total integrated bispectrum Fnl by a factor 
of 2.8 and its shape correlation with the measured bis- 



pectrum drops to 0.84, implying that the projection of 
the measured bispecrum on the tree level bispectrum is 
/ni^^ = 2.3. At larger scales the tree level bispectrum 
describes the measured bispectra much more accurately 
until late times, as expected. Similar plots for Gaussian 
A^-body simulations with 768^ and 1024^ particles are 
shown in Fig. 14d and Fig. 14e, demonstrating again the 
expected break down of tree level perturbation theory on 
small scales at late times. 



Fitting formulae for Gaussian simulations 

Separable polynomial expansion 

The separable mode expansion allowed us to com- 
press the measured A'-body dark matter bispectra to 
^max = 50 numbers f3^ at each redshift. The first ten 
of these coefficients are listed in Table II and can be used 
as a polynomial fitting formula of the matter bispectrum 
by evaluating (56). The orthogonal contributions 
appearing in (56) were defined in (51) and contain the 
nonlinear power spectrum Ps as well as the orthogonal 
polynomials Rn, which can be obtained from Qn, de- 
fined in (48), by the basis change^ Rn = ^rn^rirnQm- 



^ The transformation matrix A can be obtained from the authors 
upon request. 
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(a) Dark matter, z = A (b) Bispectrum signal, z = A 




(c) Dark matter, z = 2 (d) Bispectrum signal, z = 2 




(e) Dark matter, z = (f) Bispectrum signal, z = Q 



Figure 10. Left: Dark matter distribution in a (40Mpc//z)^ subbox of one of the C^Jq simulations at redshifts z = A^2 and 0, from top 
to bottom. Right: Measured (signal to noise weighted) bispectrum in the range 0.016/i/Mpc <k< 2/i/Mpc, averaged over the simulation 
on the left and two additional seeds. 
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Figure 11. Bispectrum slices at fixed redshift z = for different 
overall wavenumbers ksum = ki + k2 + ks. For ksum < 1/i/Mpc 
(upper row) the bispectrum signal dominates for fiattened configu- 
rations, while an enhanced equilateral contribution clearly emerges 
for ksum = 2/i/Mpc (bottom left). For ksum = 4/i/Mpc (bottom 
right) the signal is large for all triangle configurations away from 
the squeezed limit, that is, it is nearly constant. The slices show 
signal to noise weighted bispectra measured in Gaussian simula- 
tions (a) G512g and (b) 



The polynomials Qn were sorted with the slice ordering 
defined in Eq. (58) of [8]. 



Time-shift model fit 

We find a remarkably accurate fit to the matter bis- 
pectrum for Gaussian initial conditions by combining the 
modified tree level gravity shape 5f nl' defined in (8), 
with the 'constant' model co^gf defined in (14), which 
approximates the 1-halo shape, that is. 



Bf'{k,,k2,h)=Bf^l + B: 



5, const ' 



(63) 



with fitting parameters Ci and Uh- 

The combined shape (63) is dominated by the per- 
turbative gravity bispectrum at early times and by the 
constant or (approximate) halo model prediction at late 



15 



1 ^ = 10.0 


. 


11 






10 20 30 
n 

X 10^ 


40 


T 

2 = 1.0 


«n 








10 20 30 
n 

X 10^ 


40 




10 20 30 40 



500 
400 
300 
200 
100 




10 20 30 40 




10 20 30 40 
n 

(a) Gaussian 



10 20 30 40 
n 

(b) Local 



Figure 12. (a) coefficients measured with (53) for the Gaussian 
simulations G512 compared to expansion coefficients of the tree 
level gravitational bispectrum (7). (b) /S^ — (/3^) Gauss measured 
in LoclO simulations with f^f^ = 10 compared to coefficients 
of the linearly evolved local shape (see (24) and (31)). /3's were 
calculated using all modes with 0.0039/i/Mpc < k < 0.5/i/Mpc. 
For better visibility the region under the curves is colored green. 



times. We will demonstrate that this combination can 
achieve a good fit of the matter bispectrum at all red- 
shifts z < 20, while both the perturbative and the halo 
model prediction individually break down at intermedi- 
ate redshifts, when nonlinearities are important but not 
all dark matter particles can be treated as residing in 
halos. 

The fitting parameters ci and in (14) are obtained 
as follows. At each measured redshift the arbitrary 
weight w{z) in 



popt 
^6 



pgrav 



w{z){ki ^ k2 ^ ksY 



(64) 



is analytically determined such that the correlation with 
the measured bispectrum C{Bs, B^^^) is maximal (see 
green lines in Fig. 15). Then ci and are chosen such 
that CiD'^^{z) approximates the optimal weight w{z) 
over all redshifts (see black dashed lines in Fig. 15). Table 
III lists the fitting parameters obtained by this procedure. 

The shape correlation of the fitting formula (63) with 
the measured bispectrum is 99.8% or better at redshifts 
z < 20 for both /cmax = 0.5/i/Mpc and kmax = 2/i/Mpc 



22 



z = 2 



z = 



Separable expansion of theory 



Estimated from N-body 



Gravity 
(Gaussian 

initial 
conditions) 



25 30 35 



5 10 15 20 25 30 35 



Separable expansion of tlieory 



Estimated from N-body 





■ 0,4 ^-w--^, "-^ 



10 20 30 40 50 60 



10 20 30 40 50 



Separable expansion of theory 



Estimated from N-body 



Separable expansion of theory 



Estimated from N-body 



local 



floe 
NL 



10 




0.1 0.15 0.2 




05 1 15 OZ 025 3 35 



1 015 0.2 25 



equilateral 



/n1 = 100 



Separable expansion of theory Estimated from N-body 



Separable expansion of theory 



Estimated from N-body 











/I \ y. 










V 










L 


r 

01 








o.z. 







orthogonal 

/S^!' = 100 



Separable expansion of theory Estimated from N-body 




Separable expansion of theory Estimated from N-body 




flattened 



/nl = 10 



Separable expansion of theory Estimated from N-body 




Separable expansion of theory Estimated from N-body 

41 




Figure 13. Comparison of measured A/'-body bispectra (56) (on the right in each black cell) with separable expansion of tree 
level theory (44) (on the left in each black cell) for simulations G512b, LoclO, EqlOO, OrthlOO and FlatlO at redshifts z — 2 and 
z = 0. The plot axes are /ci, A:2, < 0.5/i/Mpc and we plot signal to noise weighted bispectra y^kik2k3/{Ps{ki)Ps{k2)Ps{k3))Bs 
(see Fig. 3 for comparisons at z = 30 and further descriptions of the plots). 
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Figure 14. (a)-(e): Measured bispectra in Gaussian A/'-body simulations (see Table I) as a function of redshift. Upper panels: Amplitude 
fwL^^ of the tree level gravitational bispectrum (7) (squares), normalised to unity if the tree level prediction is correct, and amplitude 
of the fitting formula from [14] (circles). Middle panels: Total integrated bispectrum (61) of measured bispectrum (symbols), tree level 
prediction (dashed lines) and fitting formula by [14] (solid lines). Lower panels: Shape correlation (59) of reconstructed bispectrum 
with separable expansions of tree level theory (squares) and fitting formula from [14] (circles). Colors indicate different k ranges for the 
bispectrum estimation. The plotted quantities are visualised in (f). 
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Figure 15. Motivation for using the growth function D in the simple fitting formula (63). 



The arbitrary weight w{^z) in S?^* 



NL +^2 + ^3)^ is determined analytically such that C{Bs, B^"^^) is maximal (for u = —1.7). We plot ||^f ^lH (black dotted), 

\\w(z)(ki + /c2 + A^s)^!! (green) and ^^st (black dashed) as defined in (14) with fitting parameters given in Table III, illustrating that 
w{z) = ciD^hi^z) is a good approximation. The continuous black and red curves show from (63) and the estimated bispectrum size 

II II, respectively. The overall normalisation can be adjusted with A^^t as explained in the main text. 
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(/3^)g512 
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FlatlO 
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^1875 
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167 
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94.16 
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1390 
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63.06 
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0.6725 


94.22 
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28.49 


362.2 


9.427 


3 


17.04 


18.58 


-183.5 


-238.7 


128.1 


0.2799 


4 


18.89 


21.7 
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-173.6 
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5 


-21.82 
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35.7 


85.29 


-43.72 


-7.911 
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9.395 


-63.53 


-6.603 


154.1 


-232.5 


1.329 
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-17.96 


-2.818 


53.24 


-77.86 


-0.4791 
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-57.18 
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171.9 


5.825 


9 


-10.25 


-21.61 


39.26 


108.1 


-109.2 


-3.548 



Table II. First 10 bispectrum expansion coefficients from 
A/'-body simulations specified in Table I for z — 0. For simu- 
lations with non-Gaussian initial conditions we show the dif- 
ference to the Gaussian simulation, (/3^)x — — (/^n )g512- 
The shape correlation between the full measured bispectrum 
(56) for 50 modes and the bispectrum corresponding to the 
shown first 10 modes is 99.7%, 99.2%, 99.95%, 95.8%, 99.7% 
and 99.9% for the columns from the left to the right. We 
used /cmax = 0.86/i/Mpc for the Gig?! simulation in the col- 
umn on the right and /cmax = 0.5/i/Mpc in all other columns. 



as shown in the upper panels of Fig. 16. The model (63), 
therefore, contains all the meaningful bispectrum shape 
information. All we require is the time dependence or 
growth rate of the bispectrum amplitude. As a first step, 
the fitting formula (64) can be normalised to the mea- 
sured bispectrum size by multiplying it with the normal- 
isation factor 
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iVfit = 



IBII 



(65) 



Table III. Fitting parameters ci and rih for the fit (63) of 
the matter bispectrum for Gaussian initial conditions (sim- 
ulations G512g and G4J0) and C2 and n^"™ for the fit (68) 
of the primordial bispectrum (57). The two columns on the 
right show the minimum shape correlation with the mea- 
sured (excess) bispectrum in A/'-body simulations, which was 
measured at redshifts z = 49, 30, 20, 10, 9, 8, . . . , 0, and the 
shape correlation at z — 0. For the equilateral case the mini- 
mum shape correlation can be improved to 99.4% if the term 
4.6 X 10-VNL5(^)°-'[2P<5(A:i)P5(fe)Fi'^(ki,k2) + 2 perms] 
is added to (68). 



which is shown by the dotted line in the lower panels of 
Fig. 16. While it varies with redshift between 0.7 and 
1.4 for /Cmax = 2/i/Mpc, it deviates by at most 8% from 
unity for /cmax = 0.5/i/Mpc. The lower panels also show 
the measured integrated bispectrum size \\B\\ and the 
two individual contributions to (63) when the normalisa- 
tion factor Nfit is included. These quantities are divided 
by II^I^j^lII for convenience. At high redshifts the total 
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(a) fcmax = 0.5/l/MpC (b) /Cmax = 2/l/MpC 

Figure 16. Upper panels: Shape correlation of estimated bispectrum in (a) G512g and (b) GHq simulations with the full 
fitting formula (63) (solid black), the shape i^f defined in (8) (dotted black) and the shape i^f ^^nst defined in (14) 

(dashed black). All correlations were calculated with 120 modes Qn for (a) k G [0.004, 0.5]/i/Mpc and (b) k G [0.016, 2]/i/Mpc. 
Lower panels: Integrated size of the measured bispectrum (red) and the contributions A^fiti^f nl (dotted black) and A^fit 5f cTnst 
(dashed black). All these bispectrum sizes are divided by the size of 5f ^1 convenience. The normalisation factor A^^t is 
defined in (65) and equals the dotted black curve. 



bispectrum size is essentially given by the contribution 
from 5f which equals the tree level prediction for the 
gravitational bispectrum in this regime. The contribu- 
tion from 5^°"^* dominates at 2; < 2 for /Cmax = 2/i/Mpc 
when filamentary and spherical nonlinear structures are 
apparent. A similar transition can be seen at later times 
on larger scales in Fig. 16a, indicating self-similar be- 
haviour. 

It is worth noting that the high integrated corre- 
lation between the simple fit (63) and measurements 
does not imply that all triangle configurations agree per- 
fectly and sub-percent level differences between shape 
correlations can in principle contain important informa- 
tion, e.g. about the observationally relevant squeezed 
limit which only makes a small contribution to the to- 
tal tetrapyd integral over the signal-to-noise weighted 
dark matter bispectrum. However, if we observed the 
dark matter bispectrum directly, these shapes would be 
hard to distinguish because the shape correlation con- 
tains the signal-to-noise weighting. Modified shape cor- 



relation weights and additional basis functions have been 
used for better quantitative comparison of the squeezed 
limit of dark matter bispectra, but this is left for a future 
publication. 

Alternative phenomenological fit 

An alternative fitting formula with 9 fitting parame- 
ters and calibrated on larger scales was given in [14] and 
summarised in this paper in Section III. In the range of 
validity given by [14], 0.03/i/Mpc < k < 0.4/i/Mpc at 
< z < 1.5, we find good agreement with our A^-body 
measurements, see green circles in Fig. 14d and Fig. 14e. 
Without having to run A'-body simulations with higher 
resolutions, we extended the bispectrum measurement to 
^max = 0.86/i/Mpc with the fast separable estimator, see 
red symbols in Fig. 14d and Fig. 14e. In this extended 
regime the fitting formula of [14] still has a shape corre- 
lation of 99.5% or more with the measured bispectrum 
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at z < 1.5, but underestimates the cumulative signal to 
noise by up to 11%. Our measured f3^ coefficients or our 
simple fit (63) can be used as alternative fitting formulae 
for the gravitational bispectrum valid to smaller scales, 
^max ^ 2/i/Mpc, and for all redshifts z <20. 

Halo model 

The halo model prediction for the Gaussian dark mat- 
ter bispectrum yields a remarkably high shape correla- 
tion of more than 99.7% with the measured bispectrum 
at z = for /Cmax = 2/i/Mpc. While the halo model bis- 
pectrum has been tested on some one- dimensional slices 
in [37] and on larger scales in [38], the result presented 
here demonstrates that at z = the halo model shape is 
a good representation of the shape measured in A/'-body 
simulations over the full tetrapyd allowed by the trian- 
gle condition. At higher redshifts, when less dark mat- 
ter resides in halos, the halo model prediction becomes 
worse and alternative phenomenological approaches like 
the ones discussed above yield higher shape correlations 
with the measured bispectrum. A more thorough exam- 
ination of the halo model will be presented in [35]. 

VII. PRIMORDIAL NON-GAUSSIAN 
BISPECTRUM RESULTS 

Primordial bispectrum measurements 

We have set up and evolved non- Gaussian initial con- 
ditions for local models with f^f^ = —10, 10, 20, 50 using 
(30), for equilateral models with = ±100, for or- 
thogonal models with /^l^ = ±100 and for the flattened 
model with /^f = 10 using separable expansions as in 
[26, 27]. Detailed parameters for the A/'-body simulations 
are given in Table I. 

Local shape 

In Fig. 12b we show comparisons of the measured 
/^n~(/^n) Gauss Coefficients with the expansion coefficients 
of the linearly evolved primordial bispectrum (24) for 
f^l = 10, finding good agreement at early times, but 
deviations at late times which are scale-dependent. The 
most obvious feature as we consider increasingly nonlin- 
ear wavenumbers is the emergence of a large constant 
or equilateral signal. Indeed Fig. 13 shows this signal on 
small scales with the late time bispectrum also exhibiting 
an enhanced signal in the squeezed limit. 

In Fig. 17a we plot the projection f^f^ of the mea- 
sured non-Gaussian bispectrum B^q on the linearly 
evolved primordial bispectrum (see (24) and (31)), their 
shape correlation C^,^, cumulative signal to noise ||5ng|| 



and total integrated bispectrum normalised to the local 
shape, ^NL*^. We also show the shape correlation Cng,g 
between ^ng and the bispectrum measured in Gaussian 
simulations. For ^ = and k G [0.0039,0.5] h/Mpc the 
linearly evolved primordial bispectrum underpredicts the 
measured cumulative signal to noise ||5ng|| by a factor 
of 8 and has a shape correlation with the measured bis- 
pectrum of about 0.6, implying that the projection f^f^ 
is about 4.8 times the input value of f^f^ = 10. The 
green and blue symbols in Fig. 17a show that tree level 
perturbation theory improves if we consider larger scales, 
as expected. 



Other shapes 

We also run simulations with equilateral, orthogonal 
and flattened non-Gaussian initial conditions using the 
separable expansion method [26, 27]. Bispectrum mea- 
surements are shown in Figures 3, 13, 17 and 18. As 
in the local case we find agreement with tree level per- 
turbation theory at all times on large scales. However 
on small scales and at late times. Fig. 13 shows that 
the measured bispectrum for equilateral initial condi- 
tions has additional equilateral and flattened contribu- 
tions compared to the linearly evolved input bispectrum 
(24). For initial conditions with the /^l^ = 100 orthog- 
onal shape, which is positive for equilateral configura- 
tions and negative for flattened and squeezed configura- 
tions, we find that the late time small scale bispectrum 
turns slightly negative for equilateral configurations and 
is more negative in flattened and squeezed configurations 
than the linearly evolved input bispectrum. The bispec- 
trum from flattened initial conditions shows additional 
contributions for squeezed, equilateral and flattened con- 
figurations compared to the tree level prediction. For 
orthogonal and flattened initial conditions the late time 
bispectra have a large signal in the squeezed limit, leading 
to possible confusion with local shape initial conditions 
(with negative f^{^ for positive input /^l^ positive 
/^IJl for positive input /nl*)- However scale dependent 
halo bias is sensitive to the scaling in the squeezed limit 
and can therefore help to disentangle such shapes (see 
e.g. [10, 24]). 

It is worth noting that the shape correlation Cf^^a be- 
tween measured and input bispectra drops as we re- 
duce /cmax already at the initial time of the simulations, 
e.g. C/3^Q, ~ 0.95 for /Cmax = 0.13/i/Mpc in the equilateral 
case. This indicates that the initial conditions generated 
with a separable mode expansion - which are generated 
at /cmax = 0.5/i/Mpc - are not a perfect fit to the shape 
when restricted to large scales. This is to be expected 
since the expansion is optimised to fit the total cumu- 
lative signal to noise up to kmax = 0.5/i/Mpc. There 
are however several possibilities of further improving the 
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Figure 17. Dark matter bispectrum quantities as in Fig. 14 but for the non-Gaussian simulations LoclO, EqlOO and OrthlOO with local, 
equilateral and orthogonal initial conditions, respectively. The non-Gaussian bispectrum is computed with (57) and then compared to the 
tree level prediction (24) (dashed lines). The lowest panel shows the shape correlation of the measured non-Gaussian bispectrum with the 
measured bispectrum for Gaussian initial conditions. 
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initial conditions on large scales. One could add basis 
functions which are localised on large scales or one could 
introduce optimised separable weighting functions in the 
mode expansions. This would be particularly important 
for examining halo bias which is very sensitive to the 
squeezed limit and therefore to large scale modes. How- 
ever for studying the global behaviour of the dark matter 
bispectrum we find it acceptable that the shape of the 
initial conditions is not perfect on very large scales and 
leave further improvements of the initial conditions for 
future work. 



on the correlation between the observations and theo- 
retical prediction Cf^^a indicate that 1-loop corrections 
can vastly improve the shape correlation in the nonlinear 
regime to 0(0.8 — 0.9), while shape correlations of more 
than 0.9 at k = 0.5/i/Mpc and z = may require higher 
order loop corrections. However we defer a detailed anal- 
ysis of loop corrections and more phenomenological halo 
model approaches to a forthcoming paper. 

Linearity in input /nl 
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Figure 18. Like Fig. 17 but for simulation Flat 10 with flat- 
tened non-Gaussian initial conditions. 



We test if corrections are important by com- 

paring simulations with f^f^ = —10,20 and 50 to the 
/nl ~ 10 simulation in Fig. 19a (with box size L = 
1600Mpc//i). The curves in all three panels would be ex- 
actly unity in case of perfect linearity in the input /nl- 
Deviations from linearity are at most 1%. 

In case of equilateral and orthogonal initial conditions 
we compare simulations with input /nl = ±100. The 
equilateral shape gives similar results to the local shape, 
see Fig. 19b. For the orthogonal shape the shape corre- 
lation between theory and measurements deviates from 
linearity in the input /nl by about 3% at z = for 
^max = 0.5/i/Mpc and the total integrated bispectrum 
differs by less than 5%. These two effects cancel each 
other approximately so that the measured projection 
/§L^ o^ly deviates by less than 1.5% from linearity. 

We conclude that at the 0(5%) level the bispectra 
measured in our large scale simulations can be scaled to 
other values f^f^ fulfilling \f^l\ < 50, \f^l\ < 100 and 
|/nl^^| ^ 100 using the linear scaling 



^ngUnl ) 



f new 
/NL 



/NL 



^ng(/nl)- 



(66) 



Fitting formulae for non- Gaussian simulations 

Separable polynomial fits 

Matter bispectra for non-Gaussian initial conditions 
of the local, equilateral, orthogonal and flattened type 
are described by the f3^ coefficients in Table II. These 
polynomial fitting formulae can serve as a starting point 
for future work that relies on non- Gaussian dark matter 
bispectra in the nonlinear regime. 



Loop corrections for primordial non-Gaussianity 

As is apparent from Figures 13, 17 and 18, the tree 
level prediction for the non- Gaussian matter bispectrum 
breaks down for small scales and low redshift. The inclu- 
sion of loop corrections described in [11] is expected to 
improve the fit at such scales and redshift s. First results 



Time-shift model fits 

Simple fitting formulae for the primordial contribution 
to the matter bispectrum can be successfully obtained 
from the halo time-shift model described in section III. 
Before discussing this fit in detail, we perform a simple 
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(a) Local (b) Equilateral (c) Orthogonal 

Figure 19. Test if the bispectrum contributions (57) from non-Gaussian initial conditions are linear in the input /nl- (a) Ratios of Loc20 
(solid line), LocSO (dash-dot line) and LoclO" (dotted line) to LoclO simulation. The label 'loc' stands for Loc20, Loc50 and LoclO". We 
scale the curves by {fl^£)in/^^ so that in case of perfect linearity in /nl all curves were unity, (b) Ratio of EqlOO" to EqlOO simulation, 
(c) Ratio of Orth~ to OrthlOO simulation. 
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Figure 20. (a) Relative growth of bispectrum size between z = 10 and ^ = as a function of the absolute value of the shape 
correlation between the measured non-Gaussian bispectrum Bng and the shape 5^°^^* defined in (14) at redshift z = 10, testing 
relation (67). The plot contains all simulations shown in (b) but data points which only differ in the input /nl are almost 
indistinguishable. Regression lines through (0, 1) are shown for /cmax = 0.5/i/Mpc (continuous) and /cmax = 2/i/Mpc (dashed), 
(b) Bispectrum size at redshift z = as a function of the bispectrum size at high redshift, z = 10, for different primordial shapes 
(see legend). Continuous lines correspond to /cmax = 0.5/i/Mpc while dashed lines correspond to /cmax = 2/i/Mpc simulations. 
Different points on one line show results for different input /nl- 



consistency check of the basic idea of the model. The 
relatively fast growth of the constant bispectrum implies 
that it constitutes the dominant contribution to the non- 
Gaussian bispectrum 5ng at sufficiently small scales and 
late time 2;iate- The amplitude of this constant bispec- 
trum is related to the projection of the non-Gaussian 
bispectrum B^q on the constant shape at the time Nearly 
when halos start to form (and thereafter). Hence we ex- 
pect that 

||^NG(^late)|| OC C(5nG (^early) , ^f, const (^early)) 

X ||5NG(^early)||. (67) 

This simple expectation is approximately seen in Fig. 20a 
for the local, equilateral and flattened shapes, confirm- 
ing the basic idea of the time-shift model. The fact that 
the orthogonal shape deviates somewhat could be related 
to the change of sign of the orthogonal shape for differ- 
ent triangle configurations. Note that relation (67) and 
Fig. 20a are interesting results on their own because they 
show that the relative growth of the non- Gaussian bis- 
pectrum can be predicted from its correlation with the 
constant shape at early times. Fig. 20b illustrates the 
absolute values of the measured bispectrum sizes which 
were used to produce Fig. 20a. 

In detail, the simple fitting formulae for the non- 
Gaussian bispectra are obtained by combining the par- 
tially loop-corrected tree level expression (25) with the 



constant shape (28) as 

fc2, fc3) ^ /nl [<T + <o"nst] • (68) 

The fitting parameters C2 and ri^^^^ in (28) are listed 
in Table IIL Similar to the Gaussian case they were 
obtained by analytically determining the optimal weight 
w{z) for the 'constant' {ki + + k^Y contribution and 
approximating this with C2D{z)^h (see green and black 
dashed lines in Fig. 21). As expected from Section III we 
find ri^^^^ = rih — I for local, equilateral and flattened 
initial conditions. 

Table III also shows the shape correlation with the 
measured excess bispectrum (57). These shape correla- 
tions are remarkably good given the simplicity of (68), 
especially for local, equilateral and flattened initial con- 
ditions. The impact of the orthogonal shape on the mat- 
ter bispectrum seems to be somewhat harder to model, 
which is reflected in the fact that our simple model per- 
forms worse for this shape. This is expected because 
the orthogonal shape contains almost no constant com- 
ponent, the dominant evolution of which is the basis 
for the simple time-shift model. While the integrated 
bispectrum size Fnl of the fitting formulae is consis- 
tent with the measured bispectrum size at high red- 
shift and at z = 0, it underestimates the measured 
size at the 10 — 20% level at intermediate redshifts for 
^max = 0.5/i/Mpc. This could be improved by adding 
more shapes or by using a redshift-dependent normalisa- 
tion factor in (68) similar to the A/fit factor in the Gaus- 
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(a) local, fcmax = 0.5/i/Mpc 



1 + z 

(b) local, fcmax = 2/i/Mpc 




^ 10' 




(c) equilateral, fcmax = 0.5/i/Mpc 



(d) equilateral, fcmax = 2/i/Mpc 



10' 





(e) flattened, fcmax = 0.5/i/Mpc 



(f) flattened, fcmax = 2/i/Mpc 




05 10' 




(g) orthogonal, fcmax = 0.5/i/Mpc 



(h) orthogonal, fcmax = 2/i/Mpc 



Figure 21. Contributions to simple non-Gaussian fits (68). The arbitrary weight w{^z) in i^^G ~ ^5 nl + ^2 + fe)^ is 

determined analytically such that C(Bng, ^ng) ^ — —1.7). We plot (black dotted) , ||il'(^)(/ci+A:2+/c3)^|| 

(green) and ^^^gt (black dashed) as defined in (28) with fitting parameters given in Table III. The continuous black and red 
curves show ||BngII from (68) and the estimated primordial bispectrum size ||Bng||, respectively. 
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sian case. On smaller scales, /Cmax = 2/i/Mpc, the non- 
Gaussian fitting formulae provide a less accurate overall 
fit, but we also list them in Table III for completeness 
and because the correlations at z = are quite high. 

It should be noted that simple fits like the ones pre- 
sented here may be somewhat more important in the 
mildly nonlinear regime than in the strongly nonlin- 
ear regime because the halo model (extended to non- 
Gaussian initial conditions) should be able to describe 
the strongly nonlinear regime. We leave more accurate 
extensions of the simple fits for non- Gaussian initial con- 
ditions and more quantitative comparisons with loop cor- 
rections and halo model predictions for future work [35] . 

Redshift at which 1-halo contribution becomes important 

We determine the redshift at which the approxi- 
mate 1-halo bispectrum 5^°"®* in the simple fits of the 
bispectrum for non-Gaussian initial conditions becomes 
important by matching 

(69) 

For /cmax = 0.5/i/Mpc we find = 3.5,3.3,4.3 for local, 
equilateral and flattened initial conditions, respectively. 
For /cmax = 2/i/Mpc we find instead z^ = 8.5,8 and 9.6 
for the same initial conditions. The orthogonal shape 
has almost no constant part. We therefore get somewhat 
different values of z^ = 5.5 and 11 for /cmax = 0.5/i/Mpc 
and /Cmax = 2/i/Mpc, respectively. This demonstrates 
that the time at which nonlinear structures contribute 
significantly to the perturbative prediction for the pri- 
mordial bispectrum is earlier on smaller scales. 

VIII. SUMMARY AND CONCLUSIONS 

We have presented an implementation of a bispec- 
trum estimator for A/'-body simulations using a separable 
modal expansion of the bispectrum as described in [26]. 
While a brute force estimation of the full bispectrum is 
computationally expensive, requiring 0{N^) operations 
for TV particles per dimension, we find that the gravita- 
tional and the most prominent primordial bispectra can 
be approximated by only nmax = 0(50) separable ba- 
sis functions (for the range of scales relevant for A/'-body 
simulations). The bispectrum projection on the corre- 
sponding subspace of all possible bispectra is estimated 
with 0(nniax^^) operations, which is faster than brute 
force estimation by a factor of 0{N^ /rima^) ~ 0(10'') 
for typical simulations. Thus the computational cost 
for accurate 3D bispectrum estimation is almost negli- 
gible compared to the cost for running the A/'-body sim- 
ulations (e.g. we can estimate the full bispectrum of a 



1024^ grid up to /Cmax = hour on only 6 

cores). This allows us to estimate the bispectrum as a 
standard simulation diagnostic whenever the power spec- 
trum is measured. Expressing the bispectrum using its 
^max separable components yields a radical compression 
of data which simplifies further analysis like comparisons 
between theory and simulations. The separable bispec- 
trum estimator therefore provides a very useful addi- 
tional statistic characterising the formation of structures 
in A/'-body simulations, with high sensitivity to different 
shapes of primordial non-Gaussianity corresponding to 
different models of inflation. 

We have performed many A/'-body simulations with 
Gaussian initial conditions as well as non-Gaussian initial 
conditions of the local, equilateral, orthogonal and (non- 
separable) flattened shape, exploiting the separable bis- 
pectrum expansion for efficient generation of initial con- 
ditions as described for primordial fields in [26, 27]. On 
large scales the measured gravitational and primordial 
bispectra agree with leading order perturbation theory 
and with measurements for Gaussian initial conditions 
by [14] in the mildly nonlinear regime, demonstrating the 
unbiasedness of the estimator and the initial conditions. 
In the nonlinear regime, the gravitational bispectrum be- 
comes dominated by a large 'constant' signal receiving 
elongated and equilateral contributions not captured by 
tree level perturbation theory. However, it remains sup- 
pressed in the squeezed limit, where primordial bispec- 
trum signals can peak (see Fig. 13). 

Our measured A/'-body bispectra for Gaussian and non- 
Gaussian simulations can be expressed by 50 components 
which we provide in Table II for the key models. 
They can be used as fitting formulae for the gravitational 
and primordial bispectrum in the nonlinear regime. 

Less accurate but simpler fitting formulae are obtained 
by modeling the bispectrum as a combination of par- 
tially loop-corrected perturbative bispectra and a sim- 
ple 'constant' bispectrum, which is constant on slices 
= const, and which is obtained as an approximation 
to the 1-halo bispectrum. While the former contribu- 
tion dominates on large scales and early times, the latter 
constant contribution dominates in the nonlinear regime. 
Interpreting the effect of primordial non-Gaussianity on 
the constant bispectrum contribution as a time-shift with 
respect to Gaussian simulations allows us to model the 
time dependence of the constant bispectrum contribu- 
tion for non- Gaussian initial conditions. For Gaussian 
initial conditions, the simple fits achieve a shape corre- 
lation of at least 99.8% with the measured gravitational 
bispectrum for z < 20 and A:niax = {0.5,2}/i/Mpc. For 
local, equilateral and flattened non-Gaussian initial con- 
ditions the primordial bispectrum is fitted with a shape 
correlation of at least 97.9% at 2; = and at least 94.4% 
for 2: < 20 (with correlation typically > 98% for most 
shapes and redshifts, see Table III). The impact of the 
orthogonal shape seems to be somewhat harder to model. 
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because it does not have a constant component initially, 
but our simple fit still achieves correlations of at least 
91% for z < 20. 

Throughout this work we have visualised the measured 
bispectra in three-dimensional tetrapyd plots [8], which 
show the bispectrum shape and amplitude at different 
length scales and generalise commonly used plots of one- 
or two-dimensional slices through the tetrapyd volume. 
For a more quantitative analysis, particularly to test an- 
alytic predictions and fitting formulae, we have made ex- 
tensive use of full three-dimensional shape correlations, 
the cumulative signal-to- noise of the bispectra and their 
projection /nl on theoretical shapes. These quantities 
have been evaluated extremely efficiently using the bis- 
pectrum components obtained from the separable esti- 
mator. 

We find that regular grid initial conditions produce an 
initial spurious bispectrum due to the anisotropy of the 
regular grid, which can be avoided by using glass initial 
conditions. However the difference between regular grid 
and glass initial conditions decreases with time as gravi- 
tational perturbations grow such that both initial condi- 
tions yield similar results at late times. Effects of order 
/nl were shown to affect the bispectrum measurements 
by less than 5% in our large scale simulations. 

Clearly, further work is required to study the effects of 
general primordial non-Gaussianity on observable quan- 
tities like the bispectrum of galaxies, particularly in the 
nonlinear regime. However the present work represents, 
we believe, an important step forward in the understand- 
ing of structure formation in the presence of primor- 
dial non-Gaussianity and the search for primordial non- 
Gaussianity in large scale structures. 
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